{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Practical Statistics for Data Scientists (Python)\n",
    "# Chapter 3. Statistial Experiments and Significance Testing\n",
    "> (c) 2019 Peter C. Bruce, Andrew Bruce, Peter Gedeck"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import required Python packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from pathlib import Path\n",
    "import random\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.stats import power\n",
    "\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "DATA = Path('.').resolve().parents[1] / 'data'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define paths to data sets. If you don't keep your data in the same directory as the code, adapt the path names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "WEB_PAGE_DATA_CSV = DATA / 'web_page_data.csv'\n",
    "FOUR_SESSIONS_CSV = DATA / 'four_sessions.csv'\n",
    "CLICK_RATE_CSV = DATA / 'click_rates.csv'\n",
    "IMANISHI_CSV = DATA / 'imanishi_data.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "session_times = pd.read_csv(WEB_PAGE_DATA_CSV)\n",
    "session_times.Time = 100 * session_times.Time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEOCAYAAABSAQgKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXAElEQVR4nO3df7DVdZ3H8edLxDQ1fwTeRTSvm9gi/UC9WzbRhLq5ZbVaY9adSsco2kkrR9fCaFdduxPOpE5u5UpLI7aiYeZoYiQZN6MdNVH8gawbKSYMipahGBnQe//4fm4ervee871wP+fn6zFzhvP9fL/fc97HL774/vx8FBGYmeWwS6MLMLP25YAxs2wcMGaWjQPGzLJxwJhZNg4YM8vGAWN1IWmlpOmNrsPqa9dGF2DtQdKmislXAy8B29L0ZyJiSv2rskaTb7Sz0SZpDfCpiPhpo2uxxvIhktWFpDWS/iG9v1DSDZL+W9ILkh6SdLik8yVtkPSkpBMq1t1H0jxJ6yWtk/RVSWMa92usLAeMNcoHgO8B+wH3Az+h+Ps4Efh34KqKZa8GtgKHAUcCJwCfqmOttoMcMNYov4iIn0TEVuAGYDwwJyK2ANcD3ZL2ldQFnAicHREvRsQG4HLgow2r3ErzSV5rlKcr3m8Gno2IbRXTAHsBBwJjgfWSBpbfBXiyHkXaznHAWLN7kuKK1Li0t2MtxIdI1tQiYj1wO3CppNdI2kXS6yW9q9G1WW0OGGsFpwG7AY8AzwE/ACY0tCIrxffBmFk23oMxs2wcMGaWjQPGzLJxwJhZNg4YM8umpW+0GzduXHR3dze6jKxefPFF9txzz0aXYTuhE7bh8uXLn42I8YPbWzpguru7uffeextdRlb9/f1Mnz690WXYTuiEbSjpiaHafYhkZtk4YMwsGweMmWXjgDGzbBwwZpaNA8bMsnHAmFk2Dhgzy6alb7Qza0YVfQcPq1P6Ycq2ByNpd0n3SHogDRt6UWq/WtLjklak19TULklXSFot6UFJR+WqzSyniNjudciXbn1FW6fIuQfzEnBcRGySNBZYJunHad55EfGDQcu/F5iUXm8Drkx/mlmLyrYHE4WB8YrHple16D4JuCatdxewryT3u2rWwrKe5JU0RtIKYAOwJCLuTrP60mHQ5ZJeldomsv1YN2tTm5m1qKwnedNAWlMl7QvcJOmNwPnAUxS9xM8FvkQxVGgpkmYCMwG6urro7+8f5aqby6ZNm9r+N3aCjt2Gg08+5XoB/wb8y6C26cCt6f1VQG/FvEeBCdU+8+ijj452tWDBgpgyZUrssssuMWXKlFiwYEGjS7IddMiXbm10CdkB98YQ/49m24ORNB7YEhF/kLQH8G7gEkkTImK9imt5JwMPp1VuAc6SdD3Fyd2NUQy61XGuu+46Zs+ezbx589i2bRtjxoxhxowZAPT29ja4OrPycp6DmQAslfQg8CuKczC3AtdKegh4CBgHfDUtfxvwGLAa+A7w2Yy1NbW+vj7mzZvHsccey6677sqxxx7LvHnz6Ovra3RpZiOSbQ8mIh4Ejhyi/bhhlg/gzFz1tJJVq1Yxbdq07dqmTZvGqlWrGlSR2Y7xowJNaPLkySxbtmy7tmXLljF58uQGVWS2YxwwTWj27NnMmDGDpUuXsnXrVpYuXcqMGTOYPXt2o0szGxE/i9SEBk7kfu5zn2PVqlVMnjyZvr4+n+C1luOAaVK9vb309vZ2RI/01r58iGRm2ThgzCwbB4yZZeOAMbNsHDBmlo0DxsyyccCYWTYOGDPLxjfaNRn3SG/txHswTWZwhz2d3CO9tT4HjJll44Axs2wcMGaWjQPGzLJxwJhZNg4YM8vGAWNm2WQLGEm7S7pH0gOSVkq6KLUfKuluSaslfV/Sbqn9VWl6dZrfnas2M6uPnHswLwHHRcRbgKnAeyQdA1wCXB4RhwHPATPS8jOA51L75Wk5M2th2QImjSi5KU2OTa8AjgN+kNrnU4zuCHBSmibNP15l7ps3s6aV9RyMpDGSVgAbgCXAb4A/RMTWtMhaYGJ6PxF4EiDN3wi8Nmd9ZpZX1ocdI2IbMFXSvsBNwN/t7GdKmgnMBOjq6qK/v39nP7LpdcJvbHedug3r8jR1RPxB0lLg7cC+knZNeykHAevSYuuAg4G1knYF9gF+N8RnzQXmAvT09ETbD+mxeJGHLWl1HbwNc15FGp/2XJC0B/BuYBWwFDglLXY6cHN6f0uaJs3/WfjRYbOWlnMPZgIwX9IYiiBbGBG3SnoEuF7SV4H7gXlp+XnA9yStBn4PfDRjbWZWB9kCJiIeBI4cov0x4K1DtP8J+HCuesys/nwnr5ll44Axs2wcMGaWjQPGzLJxwJhZNg4YM8vGAWNm2ThgzCwbB4yZZeOAMbNsHDBmlo0DxsyyccCYWTYOGDPLxgFjZtk4YMwsGweMmWXjgDGzbBwwZpZNzT55Je0OvB94J3AgsBl4GFgUESvzlmdmraxqwKQB698P9AN3U4zQuDtwODAnhc+5qYNvM7Pt1NqDuSciLhhm3mWSDgBeN9RMSQcD1wBdFGNSz42Ib0i6EPg08Exa9MsRcVta53xgBrAN+HxE/GQkP8bMmkvVgImIRYPbJO0C7BURz0fEBoq9mqFspdi7uU/S3sBySUvSvMsj4uuDPvcIirGQplAciv1U0uFp+Fkza0GlTvJKWiDpNZL2pDj/8oik86qtExHrI+K+9P4FilEdJ1ZZ5STg+oh4KSIeB1YzxPhJZtY6yl5FOiIingdOBn4MHAp8ouyXSOqmGITt7tR0lqQHJX1X0n6pbSLwZMVqa6keSGbW5MqO7DhW0liKgPlmRGyRVGrcaEl7ATcCZ0fE85KuBC6mOC9zMXAp8MmyBUuaCcwE6Orqor+/v+yqLasTfmO769RtWDZgrgLWAA8Ad0o6BHi+1koplG4Ero2IHwJExNMV878D3Jom1wEHV6x+UGrbTkTMBeYC9PT0xPTp00v+hBa1eBFt/xvbXQdvw1KHSBFxRURMjIgTo/AEcGy1dSSJYkD7VRFxWUX7hIrFPkhxTgfgFuCjkl4l6VBgEnDPCH6LmTWZWvfBnFNj/cuqzHsHxXmahyStSG1fBnolTaU4RFoDfAYgIlZKWgg8QnEF6kxfQTJrbbUOkfZOf74B+HuKvQyAD1Bj7yIilgEaYtZtVdbpA/pq1GRmLaLWfTAXAUi6EzgqXW4m3Sz3intkzMwqlb1M3QX8uWL6z6nNzGxYZa8iXQPcI+mmNH0yMD9LRWbWNkoFTET0SVoMTEtNZ0TE/fnKMrN2UHYPBmAFsH5gHUmvi4jf5ijKrJW85aLb2bh5S9VlumcNf8pynz3G8sAFJ4x2WU2hVMBI+hxwAfA0xZPOorjM/OZ8pZm1ho2bt7BmzvuGnd/f31/1Rrtq4dPqyu7BfAF4Q0T8LmcxZtZeyl5FehLYmLMQM2s/ZfdgHgP6JS0CXhporHwEwMxssLIB89v02i29zMxqKnuZeuCO3r3S9KacRZlZeyjbo90bJd0PrARWSlouaUre0sys1ZU9yTsXOCciDomIQ4Bzge/kK8vM2kHZgNkzIpYOTEREP7BnlorMrG2Uvook6V+B76Xpj1NcWTIzG1bZPZhPAuOBH1J0gTmOEfSja2adqexVpOeAz2euxczaTNmrSEsk7VsxvZ8kj7poZlWVPUQaFxF/GJhIezQHZKnIzNpG2YD5i6S/jkGdhi0pNS6SmXWusleRZgPLJP2coquGd5IGPzMzG07ZcZEWA0cB3weuB46OiKrnYCQdLGmppEckrZT0hdS+fzqn8+v0536pXZKukLQ6DSt71M79NDNrtLIneQW8h2JkgVuBV0uqNTD9VuDciDgCOAY4U9IRwCzgjoiYBNyRpgHeSzHY2iSKvaMrR/pjzKy5lD0H823g7UBvmn4B+Fa1FSJifUTcl96/AKyiGMz+JF7uMHw+RQfipPZr0siRdwH7DhoF0sxaTNmAeVtEnAn8Cf56Fal0tw2SuoEjgbuBrohYn2Y9xcvDn0yk6NhqwNrUZmYtquxJ3i2SxpCuHEkaD/ylzIqpi4cbgbMj4vniaKsQESFpRFejJM0knWDu6uqiv79/JKu3pE74ja2u2jbatGlTzW3Yrtu4bMBcAdwEHCCpDzgF+EqtlSSNpQiXayPih6n5aUkTImJ9OgTakNrXAQdXrH5QattORMyleLqbnp6eqNaZcltYvKhqh9HWBGpso1qdfrfzNi57Fela4IvA1yiGLjk5Im6otk46MTwPWDWoa81bgNPT+9OBmyvaT0tXk44BNlYcSplZCyo7bMnrgccj4luSpgPvlrS+8u7eIbwD+ATwkKQVqe3LwBxgoaQZwBPAqWnebcCJwGrgj8AZI/olZtZ0yh4i3Qj0SDoMuIpib2MBRSAMKSKWUdyUN5Tjh1g+gDNL1mNmLaD0owIRsRX4EPDNiDgP8CVkM6uqbMBskdQLnAbcmtrG5inJzNpF2YA5g+JGu76IeFzSobzcu52Z2ZDKdjj1CBUdTkXE48AluYoys/ZQdg/GzGzEHDBmlo0DxsyyKXuj3eHAecAhletExHGZ6jKzNlD2RrsbgP+kGM1xW75yzKydlA2YrRHhDqDMbETKnoP5kaTPSpqQurzcX9L+WSszs5ZXdg9m4Onn8yraAvjb0S2ns7zlotvZuHlLzeW6Zy0adt4+e4zlgQtOGM2yzEZN2RvtDs1dSCfauHkLa+a8r+oytfoSqRY+Zo1WNWAkHRcRP5P0oaHmV3QiZWb2CrX2YN4F/Az4wBDzAnDAmNmwqgZMRFyQ/nTnT2Y2YlWvIkn6uKRhl5H0eknTRr8sM2sHtQ6RXgvcL2k5sBx4BtgdOIzi8OlZXh44zcxsO7UOkb4h6ZvAcRR97L4Z2EwxiNonIuK3+Us0s1ZV8zJ1RGwDlqSXmVlpfprazLLJFjCSvitpg6SHK9oulLRO0or0OrFi3vmSVkt6VNI/5qrLzOon5x7M1cB7hmi/PCKmptdtAJKOAD4KTEnrfDsNVWtmLaxUwEjqkjRP0o/T9BFp4LRhRcSdwO9L1nEScH1EvJT6+10NvLXkumbWpMruwVwN/AQ4ME3/H3D2Dn7nWZIeTIdQ+6W2icCTFcusTW1m1sLKPk09LiIWSjofICK2StqRjqeuBC6meMzgYuBS4JMj+QBJM4GZAF1dXfT39+9AGc2jVv2bNm2quUyr/zdoB9W2QSdvw7IB86Kk11IEAwOD04/0yyLi6YH3kr7Dy4O4rQMOrlj0oNQ21GfMBeYC9PT0RLUnjZve4kVVn5SG2k9Tl/kMy6zGNujkbVg2YM6hGI/69ZJ+CYwHThnpl0maEBHr0+QHgYErTLcACyRdRnEYNgm4Z6Sfb9YIe0+exZvm17ihfX619QGqd9vRqsr2B3OfpHcBb6AY0P7RiKjaU5Kk64DpwDhJa4ELgOmSplLsCa0BPpM+f6WkhcAjwFbgzHSDn1nTe2HVnKr9+nRynz5lRxUYA5wIdKd1TpBERFw23DoR0TtE87wqy/cBfWXqMbPWUPYQ6UfAn4CHgL/kK8fM2knZgDkoIt6ctRIzaztl74P5sST3LG1mI1J2D+Yu4KbU+dQWihO9ERGvyVaZmbW8sgFzGfB24KGIiIz1mFkbKXuI9CTwsMPFzEai7B7MY0B/etjxpYHGapepzczKBszj6bVbepmZ1VT2Tt6LchdiZu2n1siO34yIsyT9iPSgY6WI+KdslZlZy6u1B3MacBbw9TrUYmZtplbA/AYgIn5eh1rMrM3UCpjxks4ZbqavIplZNbUCZgywF8Wdu2ZmI1IrYNZHxL/XpRIzazu17uT1nouZ7bBaAXN8Xaows7ZUNWAiouy4RmZmr+Cxqc0sGweMmWXjgDGzbLIFTBoadoOkhyva9pe0RNKv05/7pXZJukLS6jSs7FG56jKz+sm5B3M18J5BbbOAOyJiEnBHmgZ4L8Vga5MohoW9MmNdZlYn2QImIu4EBl+FOomXx7ibD5xc0X5NFO4C9pU0IVdtZlYf9T4H01UxdOxTQFd6P5GiW84Ba1ObmbWwsj3ajbqICEkj7uNX0kyKwyi6urro7+8f7dLqqlb9mzZtqrlMq/83aAfVtkEnb8N6B8zTkiZExPp0CLQhta8DDq5Y7qDU9goRMReYC9DT0xPVxvxteosXVR2zGGqPa1zmMyyzGtugk7dhvQ+RbgFOT+9PB26uaD8tXU06BthYcShlZi0q2x6MpOuA6cA4SWuBC4A5wEJJM4AngFPT4rcBJwKrgT8CZ+Sqq5nsPXkWb5o/q/aC84eftfdkgPeNVklmoypbwERE7zCzXvEAZRpv6cxctTSrF1bNYc2c6uFQa/e6e9aiUa7KbPT4Tl4zy8YBY2bZOGDMLBsHjJll44Axs2wcMGaWjQPGzLJp2LNIZu2k5v1Ii4efv88eY0e5mubhgDHbSbVuluyetajmMu3Kh0hmlo0DxsyyccCYWTYOGDPLxgFjZtk4YMwsG1+mbrBS/bl06D0U1vocMA1U5t6ITr6HwlqfD5HMLBsHjJll44Axs2wcMGaWTUNO8kpaA7wAbAO2RkSPpP2B7wPdwBrg1Ih4rhH1mdnoaOQezLERMTUietL0LOCOiJgE3JGmzayFNdMh0km8PMTYfODkxpViZqOhUQETwO2SlqfB7AG6KoaLfQroakxpZjZaGnWj3bSIWCfpAGCJpP+tnBkRISmGWjEF0kyArq4u+vv7sxfbaJ3wG9tdp27DhgRMRKxLf26QdBPwVuBpSRMiYr2kCcCGYdadC8wF6OnpiWrDqraFxYuqDh1rLaCDt2HdD5Ek7Slp74H3wAnAw8AtwOlpsdOBm+tdm5mNrkbswXQBN0ka+P4FEbFY0q+AhZJmAE8ApzagNjMbRXUPmIh4DHjLEO2/A46vdz1moy3947l92yXbT0cMeYqx7TTTZWqzthAR272WLl36irZO4YAxs2wcMGaWjQPGzLJxwJhZNg4YM8vGAWNm2ThgzCwbB4yZZeOAMbNsHDBmlo0DxsyyccCYWTYeOrbJ+Elcayfeg2kyfhLX2okDxsyyccCYWTYOGDPLxgFjZtk4YMwsGweMmWXjgDGzbBwwZpaNWvnGLUnPUAzS1s7GAc82ugjbKZ2wDQ+JiPGDG1s6YDqBpHsjoqfRddiO6+Rt6EMkM8vGAWNm2Thgmt/cRhdgO61jt6HPwZhZNt6DMbNsHDCZSdomaYWkhyXdIOnVmb9vnKQtkv455/d0knpuQ0n9kh5N37dK0sxc31UPDpj8NkfE1Ih4I/BnIPf/+B8G7gJ6M39PJ6n3NvxYREwF3gFcImm3zN+XjQOmvn4BHCbpA5LulnS/pJ9K6gKQNF7SEkkrJf2XpCckjUvzPi7pnvQv21WSxgzzHb3AucBESQfV52d1lHpswwF7AS8C2/L+pHwcMHUiaVfgvcBDwDLgmIg4Erge+GJa7ALgZxExBfgB8Lq07mTgI8A70r9s24CPDfEdBwMTIuIeYGFax0ZJPbZhcq2kB4FHgYsjomUDxp1+57eHpBXp/S+AecAbgO9LmgDsBjye5k8DPggQEYslPZfajweOBn6VOgXfA9gwxHd9hCJYoPhL/13g0tH8MR2qntsQikOkeyWNB/5H0uKIaMlHYhww+W1O/2L9laT/AC6LiFskTQcurPEZAuZHxPk1lusF/kbSwL+MB0qaFBG/HnHVVqme2/CvIuIZSfcBb6NFn7nzIVJj7AOsS+9Pr2j/JXAqgKQTgP1S+x3AKZIOSPP2l3RI5QdKOhzYKyImRkR3RHQDX8Mne3MZ9W04WLpadSTwm1Gsu64cMI1xIXCDpOVs/5TtRcAJkh6muBr0FPBCRDwCfAW4PR2bLwEmDPrMXuCmQW034oDJ5UJGfxsOuDYdki0Hro6I5Xl+Qn6+k7eJSHoVsC0itkp6O3Dl4F1za27ehtvzOZjm8jpgoaRdKO63+HSD67GR8zas4D0YM8vG52DMLBsHjJll44Axs2wcMGaWjQPGzLJxwJhZNv8Po6FuQ3ux/poAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = session_times.boxplot(by='Page', column='Time',\n",
    "                           figsize=(4, 4))\n",
    "ax.set_xlabel('')\n",
    "ax.set_ylabel('Time (in seconds)')\n",
    "plt.suptitle('')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "35.66666666666667\n"
     ]
    }
   ],
   "source": [
    "mean_a = session_times[session_times.Page == 'Page A'].Time.mean()\n",
    "mean_b = session_times[session_times.Page == 'Page B'].Time.mean()\n",
    "print(mean_b - mean_a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code is different to the R version. idx_A and idx_B are reversed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "61.266666666666666\n"
     ]
    }
   ],
   "source": [
    "# Permutation test example with stickiness\n",
    "def perm_fun(x, nA, nB):\n",
    "    n = nA + nB\n",
    "    idx_B = set(random.sample(range(n), nB))\n",
    "    idx_A = set(range(n)) - idx_B\n",
    "    return x.loc[idx_B].mean() - x.loc[idx_A].mean()\n",
    "    \n",
    "nA = session_times[session_times.Page == 'Page A'].shape[0]\n",
    "nB = session_times[session_times.Page == 'Page B'].shape[0]\n",
    "print(perm_fun(session_times.Time, nA, nB))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgX0lEQVR4nO3de7xVZb3v8c9XUC6amoFuMhP1oG0IIVhYbW+oGSrmXcm8V9syLcy2Rwwya2svu3rEY5mpW0kSxBJN7IL3C6mxCLmoeElQlJSkFJQjIL/zx3gWTldrLebCNeazWOv7fr3mizGeOeYYvznmXF/GfOYYz1REYGZmtbdJ7gLMzDorB7CZWSYOYDOzTBzAZmaZOIDNzDLpmruA96JXr17Rt2/f3GWYZVdfXw/A0KFDM1diTamvr/97RPRu3K6N+TS0urq6mDlzZu4yzLKTBMDG/PfckUmqj4i6xu3ugjAzy8QBbGaWiQPYzCwTB7CZWSYOYDOzTBzAZmaZOIDNzDJxAJuZZeIANjPLxAFsZpaJA9jMLBMHsJlZJg5gM7NMNurhKM2q0XfMtNLWvfCSkaWt2zo+HwGbmWXiADYzy8QBbGaWiQPYzCwTB7CZWSYOYDOzTBzAZmaZOIDNzDJxAJuZZeIANjPLxAFsZpaJA9jMLBMHsJlZJg5gM7NMHMBmZpk4gM3MMnEAm5ll4gA2M8vEAWxmlokD2MwsEwewmVkmDmAzs0wcwGZmmTiAzQyAvn37IqlT3Pr27Zt7dwPQNXcBZtY+LFq0iIjIXUZNSMpdAuAjYDOzbBzAZmaZOIDNzDJxAJtZsxYvXszhhx9Ov3792GWXXRg9ejSrVq3iuuuu46yzzspd3r/YYostcpfQKg5gM2tSRHDUUUdxxBFH8PTTT/PUU0+xYsUKxo4dW8r21qxZU8p62zMHsJk16e6776Z79+6cdtppAHTp0oVLL72Ua6+9ljfffJMXXniB4cOH069fP77zne8A8MYbbzBy5EgGDRrERz/6USZPngxAfX09++67L0OHDmXEiBEsWbIEgOHDh3P22WdTV1fHxRdfzI477sjatWvXrWuHHXZg9erVPPvssxx00EEMHTqUvffemyeffBKA5557jk9+8pMMHDiQcePG1XoXvWelnYYmaQdgArAdEMBVEXGZpG2AyUBfYCFwXET8Q8V5IZcBhwBvAqdGxKyy6rP2oe+YaaWte+ElI0tbd2cwf/58hg4d+q62Lbfckg9/+MOsWbOGRx99lHnz5tGzZ0+GDRvGyJEjWbRoER/84AeZNq14XV977TVWr17NV7/6VW699VZ69+7N5MmTGTt2LNdeey0Aq1atYubMmQDMmjWL++67j/3224/bb7+dESNGsOmmm3L66adz5ZVX0q9fPx555BG+8pWvcPfddzN69GjOOOMMTj75ZK644ora7qA2UOYR8BrgGxHRH/gEcKak/sAY4K6I6AfcleYBDgb6pdvpwM9KrM3M3qMDDzyQD3zgA/To0YOjjjqKBx98kIEDBzJ9+nTOO+88HnjgAbbaaisWLFjAvHnzOPDAAxk8eDAXXXQRixcvXreeUaNGvWu64ah50qRJjBo1ihUrVjBjxgyOPfZYBg8ezJe+9KV1R9APPfQQxx9/PAAnnXRSDZ992yjtCDgilgBL0vRySU8A2wOHA8PTYtcD9wLnpfYJUZwJ/rCkrSX1Sesxsxrr378/N99887vaXn/9dZ5//nm6du36LxczSGLXXXdl1qxZ3HHHHYwbN44DDjiAI488kgEDBvCnP/2pye1svvnm66YPO+wwvvnNb7Js2TLq6+vZf//9eeONN9h6662ZPXt2k49vLxdVbIia9AFL6gt8DHgE2K4iVP9G0UUBRTi/UPGwxamt8bpOlzRT0sylS5eWV7RZJ3fAAQfw5ptvMmHCBADefvttvvGNb3DqqafSs2dPpk+fzrJly1i5ciVTp05lzz335KWXXqJnz56ceOKJnHvuucyaNYvddtuNpUuXrgvg1atXM3/+/Ca3ucUWWzBs2DBGjx7NoYceSpcuXdhyyy3ZaaedmDJlClB8OfjYY48BsOeeezJp0iQAJk6cWPYuaXOlB7CkLYBfA2dHxOuV96Wj3VZd+xgRV0VEXUTU9e7duw0rNbNKkrjllluYMmUK/fr1Y9ddd6V79+5873vfA2CPPfbg6KOPZvfdd+foo4+mrq6OuXPnssceezB48GC+853vMG7cODbbbDNuvvlmzjvvPAYNGsTgwYOZMWNGs9sdNWoUN9xww7u6JiZOnMg111zDoEGDGDBgALfeeisAl112GVdccQUDBw7kxRdfLHeHlEBlXvstaVPgduAPEfGT1LYAGB4RSyT1Ae6NiN0k/TxN39h4uebWX1dXFw2d97ZxqsWXcJ3hi76Gj+Hv5e9ZUqcaC6KWz1VSfUTUNW4v7Qg4ndVwDfBEQ/gmtwGnpOlTgFsr2k9W4RPAa+7/NbOOrMzR0PYETgLmSpqd2r4JXALcJOkLwCLguHTfHRSnoD1DcRraaSXWZmaWXZlnQTwINPf15AFNLB/AmWXVY2bW3vhKODOryoUXXsiPfvQjLrjgAu68804AHnjgAQYMGMDgwYNZuXIl5557LgMGDODcc8/NXO3GwQOym1mrfPe73103PXHiRM4//3xOPPFEAK666iqWLVtGly5dqlrXmjVr6Nq188aQj4DNrFkXX3wxu+66K3vttRcLFiwA4NRTT+Xmm2/m6quv5qabbuJb3/oWJ5xwAocddhgrVqxg6NChTJ48maVLl3L00UczbNgwhg0bxkMPPQQUR9InnXQSe+65JyeddFKLy33+859n+PDh7LzzzowfP35dXRMmTGD33Xdn0KBB666Aa2497Vnn/a/HzFpUX1/PpEmTmD17NmvWrGHIkCHvGhvii1/8Ig8++CCHHnooxxxzDFBcSNFwxdrnPvc5vv71r7PXXnvx/PPPM2LECJ544gkAHn/8cR588EF69OjR4nJPPvkk99xzD8uXL2e33XbjjDPO4KmnnuKiiy5ixowZ9OrVi2XLlgEwevToZtfTXjmAzaxJDzzwAEceeSQ9e/YEisuEW+POO+/k8ccfXzf/+uuvs2LFinXr6tGjx3qXGzlyJN26daNbt25su+22vPzyy9x9990ce+yx9OrVC4BtttmmxfW05zGCHcBmVoq1a9fy8MMP071793+5r3L8h5aW69at27rpLl26tDhmcEvraa/cB2xmTdpnn32YOnUqK1euZPny5fz2t79t1eM//elPc/nll6+bb24wnWqXa7D//vszZcoUXn31VYB1XRCtXU974AA2syYNGTKEUaNGMWjQIA4++GCGDRvWqsePHz+emTNnsvvuu9O/f3+uvPLK97RcgwEDBjB27Fj23XdfBg0axDnnnLNB62kPSh0LomweC2Lj57Eg2obHgmidDj8WhJmZtcwBbGaWiQPYzCwTB7CZWSY+D9jMANhxxx036t9Xa40dd9wxdwmAA9jMkoULF+YuodNxF4SZWSYOYDOzTBzAZmaZOIDNzDJxAJuZZeIANjPLxAFsZpaJA9jMLBMHsJlZJg5gM7NMHMBmZpk4gM3MMvFgPGZtoKyfPWovP3lk5fARsJlZJg5gM7NMHMBmZpk4gM3MMnEAm5ll4gA2M8vEAWxmlokD2MwsEwewmVkmDmAzs0wcwGZmmTiAzcwycQCbmWXiADYzy8QBbGaWiQPYzCwTB7CZWSYOYDOzTBzAZmaZOIDNzDJxAJuZZeIANjPLxAFsZpaJA9jMLBMHsJlZJg5gM7NMHMBmZpmUFsCSrpX0iqR5FW0XSnpR0ux0O6TivvMlPSNpgaQRZdVlZtZelHkEfB1wUBPtl0bE4HS7A0BSf+CzwID0mJ9K6lJibWZm2ZUWwBFxP7CsysUPByZFxFsR8RzwDLBHWbWZmbUHOfqAz5I0J3VRvD+1bQ+8ULHM4tT2LySdLmmmpJlLly4tu1Yzs9LUOoB/BuwCDAaWAD9u7Qoi4qqIqIuIut69e7dxeWZmtVPTAI6IlyPi7YhYC/yCd7oZXgR2qFj0Q6nNzKzDqmkAS+pTMXsk0HCGxG3AZyV1k7QT0A94tJa1mZnVWteyVizpRmA40EvSYuDbwHBJg4EAFgJfAoiI+ZJuAh4H1gBnRsTbZdVmZtYelBbAEXF8E83XtLD8xcDFZdVjZtbe+Eo4M7NMHMBmZpk4gM3MMnEAm5ll4gA2M8vEAWxmlokD2MwsEwewmVkmDmAzs0wcwGZmmTiAzcwyqSqAJQ0suxAzs86m2iPgn0p6VNJXJG1VakVmZp1EVQEcEXsDJ1AMml4v6VeSDiy1MjOzDq7qPuCIeBoYB5wH7AuMl/SkpKPKKs7MrCOrtg94d0mXAk8A+wOfiYh/T9OXllifmVmHVe2A7JcDVwPfjIiVDY0R8ZKkcaVUZmbWwVUbwCOBlQ0/EyRpE6B7RLwZEb8srTozsw6s2j7gO4EeFfM9U5uZmW2gagO4e0SsaJhJ0z3LKcnMrHOoNoDfkDSkYUbSUGBlC8ubmdl6VNsHfDYwRdJLgIB/A0aVVZS1D33HTCtt3QsvGVnaus02FlUFcET8WdJHgN1S04KIWF1eWWZmHV+1R8AAw4C+6TFDJBERE0qpysysE6gqgCX9EtgFmA28nZoDcACbmW2gao+A64D+ERFlFmNm1plUexbEPIov3szMrI1UewTcC3hc0qPAWw2NEXFYKVWZmXUC1QbwhWUWYWbWGVV7Gtp9knYE+kXEnZJ6Al3KLc3MrGOrdjjK/wRuBn6emrYHppZUk5lZp1Dtl3BnAnsCr8O6wdm3LasoM7POoNoAfisiVjXMSOpKcR6wmZltoGoD+D5J3wR6pN+CmwL8tryyzMw6vmoDeAywFJgLfAm4g+L34czMbANVexbEWuAX6WZmZm2g2rEgnqOJPt+I2LnNKzIz6yRaMxZEg+7AscA2bV+OmVnnUVUfcES8WnF7MSL+D8UPdZqZ2QaqtgtiSMXsJhRHxK0ZS9jMzBqpNkR/XDG9BlgIHNfm1ZiZdSLVngWxX9mFmJl1NtV2QZzT0v0R8ZO2KcfMrPNozVkQw4Db0vxngEeBp8soysysM6g2gD8EDImI5QCSLgSmRcSJZRVmZtbRVXsp8nbAqor5VanNzMw2ULVHwBOARyXdkuaPAK4vpSIzs06i2rMgLpb0O2Dv1HRaRPylvLLMzDq+arsgAHoCr0fEZcBiSTuVVJOZWadQ7U8SfRs4Dzg/NW0K3FBWUWZmnUG1R8BHAocBbwBExEvA+8oqysysM6g2gFdFRJCGpJS0eXklmZl1DtWeBXGTpJ8DW6dfSP48HpzdrKb6jpnWJss0tvASD2yYy3oDWJKAycBHKH4VeTfggoiYXnJtZmYd2noDOCJC0h0RMRBw6JqZtZFq+4BnSRpWaiVmZp1MtQH8ceBhSc9KmiNprqQ5LT1A0rWSXpE0r6JtG0nTJT2d/n1/apek8ZKeSesf0vyazcw6hhYDWNKH0+QIYGdgf4qR0A5N/7bkOuCgRm1jgLsioh9wV5oHOBjol26nAz+rrnwzs43X+o6ApwJExCLgJxGxqPLW0gMj4n5gWaPmw3lnDInrKcaUaGifEIWHKc626FP90zAz2/isL4BVMd0WP0G/XUQsSdN/450R1bYHXqhYbnFq+9eCpNMlzZQ0c+nSpW1QkplZHusL4Ghm+j2rvLCjlY+7KiLqIqKud+/ebVmSmVlNre80tEGSXqc4Eu6RpknzERFbtnJ7L0vqExFLUhfDK6n9RWCHiuU+lNrMzDqsFo+AI6JLRGwZEe+LiK5pumG+teELxU8anZKmTwFurWg/OZ0N8QngtYquCjOzDqnaS5FbTdKNwHCgl6TFwLeBSygua/4CsIh3ftr+DuAQ4BngTeC0suoyM2svSgvgiDi+mbsOaGLZAM4sqxYzs/aoNQOym5lZG3IAm5ll4gA2M8vEAWxmlokD2MwsEwewmVkmDmAzs0wcwGZmmTiAzcwycQCbmWXiADYzy8QBbGaWiQPYzCwTB7CZWSYOYDOzTBzAZmaZOIDNzDJxAJuZZeIANjPLxAFsZpaJA9jMLBMHsJlZJg5gM7NMHMBmZpk4gM3MMnEAm5ll4gA2M8vEAWxmlokD2MwsEwewmVkmDmAzs0wcwGZmmTiAzcwycQCbmWXiADYzy8QBbGaWiQPYzCwTB7CZWSYOYDOzTBzAZmaZOIDNzDJxAJuZZeIANjPLxAFsZpaJA9jMLJOuuQuwDdd3zLRS1rvwkpGlrNfM3s1HwGZmmTiAzcwycQCbmWXiADYzy8QBbGaWiQPYzCwTB7CZWSYOYDOzTBzAZmaZZLkSTtJCYDnwNrAmIuokbQNMBvoCC4HjIuIfOeozM6uFnEfA+0XE4IioS/NjgLsioh9wV5o3M+uw2lMXxOHA9Wn6euCIfKWYmZUvVwAH8EdJ9ZJOT23bRcSSNP03YLumHijpdEkzJc1cunRpLWo1MytFrtHQ9oqIFyVtC0yX9GTlnRERkqKpB0bEVcBVAHV1dU0uY2a2MchyBBwRL6Z/XwFuAfYAXpbUByD9+0qO2szMaqXmR8CSNgc2iYjlafrTwHeB24BTgEvSv7fWujazzs5jTNdWji6I7YBbJDVs/1cR8XtJfwZukvQFYBFwXIbazMxqpuYBHBF/BQY10f4qcECt6zEzy6U9nYZmZtapOIDNzDJxAJuZZeIANjPLxAFsZpaJA9jMLBMHsJlZJg5gM7NMHMBmZpk4gM3MMnEAm5ll4gA2M8vEAWxmlokD2MwsEwewmVkmDmAzs0wcwGZmmTiAzcwycQCbmWXiADYzy8QBbGaWiQPYzCwTB7CZWSZdcxfQUfUdM62U9S68ZGQp6zWz2vMRsJlZJg5gM7NMHMBmZpk4gM3MMnEAm5ll4gA2M8vEAWxmlokD2MwsEwewmVkmDmAzs0wcwGZmmTiAzcwycQCbmWXiADYzy8QBbGaWiccDNrOa8ljZ7/ARsJlZJg5gM7NMHMBmZpk4gM3MMnEAm5ll4gA2M8vEAWxmlokD2MwsEwewmVkmDmAzs0wcwGZmmXTKsSB8LbpZx7ax/I37CNjMLBMHsJlZJg5gM7NMHMBmZpm0uwCWdJCkBZKekTQmdz1mZmVpVwEsqQtwBXAw0B84XlL/vFWZmZWjXQUwsAfwTET8NSJWAZOAwzPXZGZWCkVE7hrWkXQMcFBEfDHNnwR8PCLOqljmdOD0NLsbsKCJVfUC/l5yudVyLU1zLU1zLU3b2GvZMSJ6N27c6C7EiIirgKtaWkbSzIioq1FJLXItTXMtTXMtTeuotbS3LogXgR0q5j+U2szMOpz2FsB/BvpJ2knSZsBngdsy12RmVop21QUREWsknQX8AegCXBsR8zdgVS12UdSYa2maa2maa2lah6ylXX0JZ2bWmbS3Lggzs07DAWxmlslGH8CSjpU0X9JaSXWN7js/XdK8QNKIivbSL3eWNFnS7HRbKGl2au8raWXFfVeWsf1GtVwo6cWKbR5ScV+T+6jEWn4o6UlJcyTdImnr1F7z/ZK2m+XSd0k7SLpH0uPp/Ts6tTf7WpVcz0JJc9M2Z6a2bSRNl/R0+vf9Nahjt4rnPlvS65LOrtV+kXStpFckzatoa3I/qDA+vXfmSBrS6g1GxEZ9A/6d4oKMe4G6ivb+wGNAN2An4FmKL/a6pOmdgc3SMv1LrvHHwAVpui8wr8b76ELgv5pob3IflVzLp4Guafr7wPcz7peavxcqtt0HGJKm3wc8lV6PJl+rGtSzEOjVqO0HwJg0Pabhtarx6/M3YMda7RdgH2BI5Xuxuf0AHAL8DhDwCeCR1m5voz8CjognIqKpq+EOByZFxFsR8RzwDMWlzjW93FmSgOOAG8vaxnvQ3D4qTUT8MSLWpNmHKc71ziXbpe8RsSQiZqXp5cATwPa12HYrHA5cn6avB46o8fYPAJ6NiEW12mBE3A8sa9Tc3H44HJgQhYeBrSX1ac32NvoAbsH2wAsV84tTW3PtZdkbeDkinq5o20nSXyTdJ2nvErdd6az0Menaio+Std4XjX2e4giiQa33S+7nDxTdL8DHgEdSU1OvVdkC+KOk+nS5P8B2EbEkTf8N2K5GtTT4LO8+cMmxX6D5/fCe3z8bRQBLulPSvCZuWQfqqbKu43n3m2gJ8OGI+BhwDvArSVuWXMvPgF2AwWn7P36v23sPtTQsMxZYA0xMTaXsl/ZO0hbAr4GzI+J1avxaVdgrIoZQjER4pqR9Ku+M4jN3zc5ZVXEh1mHAlNSUa7+8S1vvh3Z1IUZzIuJTG/Cwli5rbpPLnddXl6SuwFHA0IrHvAW8labrJT0L7ArM3JAaqq2loqZfALen2VIu/a5iv5wKHAockN7Qpe2X9ch66bukTSnCd2JE/AYgIl6uuL/ytSpVRLyY/n1F0i0U3TMvS+oTEUvSR+tXalFLcjAwq2F/5NovSXP74T2/fzaKI+ANdBvwWUndJO0E9AMepbaXO38KeDIiFjc0SOqtYtxjJO2c6vprSdtv2GZlv9SRQMM3vM3tozJrOQj438BhEfFmRXvN9wsZL31P3w1cAzwRET+paG/utSqzls0lva9hmuKL0nkU++KUtNgpwK1l11LhXZ8cc+yXCs3th9uAk9PZEJ8AXqvoqqhOLb/VLOlbyyMp+l7eAl4G/lBx31iKb7kXAAdXtB9C8a3zs8DYEmu7Dvhyo7ajgfnAbGAW8Jka7KNfAnOBOelN02d9+6jEWp6h6DebnW5X5tovtXwvNLHdvSg+ys6p2BeHtPRalVjLzhRngDyWXoOxqf0DwF3A08CdwDY12jebA68CW1W01WS/UIT+EmB1ypUvNLcfKM5+uCK9d+ZScRZWtTdfimxmlklH7oIwM2vXHMBmZpk4gM3MMnEAm5ll4gA2M8vEAdzOSRqrYrSsOWkUqI+30Xo/KOnmNlrX2ZJ6VszfoTTKWVkkXafiV7SRdLWk/mn6WElPSLonzd+Y9t3Xy6ynrUnqI+n2NF0naXzumlpDxeh2zZ6rK2kzSfeni5U6rU795Ns7SZ+kuGJsSES8JakXxahd71lEvAQc0xbrAs4GbgDeTOuuyRCKDSLiixWzXwD+MyIelPRvwLCI+F/VrktS13hnsKCczgF+ARARMyn3isCai4hVku4CRvHO5eidjo+A27c+wN+juEyXiPh7Ck4kDU2D1tRL+kPDlUKSvqZijNk5kialtn31zjiqf5H0vsojFEndJf2PivFg/yJpv9R+qqTfSPq9irFQf9C4QElfAz4I3FNx1LlQUq+0jSfT0epTkiZK+pSkh9L69kjLb65igJVH0/b/ZYyPdLXR/1Uxdu+dwLYV992bjhIvoLjA4RpJPwT+CGyfnvfeknZJz6Ve0gOSPpIef52kKyU9AvxgPcuNlzRD0l8bjsDTfeel/feYpEtSW3PrOVbF2BiPSbq/mdf+aOD3afnhFUfDF6Z9dW+q4WtN7KsuqdZ5qaavr6ee7VSMzfxYuv1Haj9H74zjcXZq66viE8YvVHwy+6OkHhXvycckPQacWVHPgPTazk7vy37prqnACc08/86hVlf++LZBV+VsQXGF1FPAT4F9U/umwAygd5ofRfEDpgAvAd3S9Nbp398Ce1assysV4+8C36h4/EeA54HuwKkUlwNvleYXATs0UedCKsaSbZhP21gDDKT4z74euJbiCqLDgalp+e8BJzbUnJ7v5o22cRQwnWKM2A8C/wSOSffdS7oKqdH0uueY5u8C+qXpjwN3p+nrKMYW6FLFclPSc+lPMZQlFOMWzAB6pvlt1rOeucD2la9Ro+e6E1BfMT8cuD1NX5i21S3t41eBTRs9figwvWJ+6/XUM5liMCDS/t0qrWMuxVVpW1BcIfexitd0cFr+porXbg6wT5r+Ie+8vy4HTkjTmwE9Kra1NPffWc6buyDasYhYIWkoxZCW+wGTVfxqw0zgo8B0SVC8kRuuQZ8DTJQ0leIIA+Ah4CeSJgK/iYjF6XEN9qL4IyEinpS0iGIgHIC7IuI1AEmPUwyOXTkE3/o8FxFz0+Pnp/WFpLkUf8xQjD1wmKT/SvPdgQ9TjJHbYB/gxoh4G3hJ0t2tqKFh1LH/AKZUPPduFYtMiYi3q1huakSsBR6X1DAs4aeA/4k0tkVELFvPeh4CrpN0E/CbJsrtAyxt4elMizR4kaRXKIZHXFxx/1+BnSVdDkyjGGaypXr2B05Otb8NvCZpL+CWiHgDQNJvKN6Ht1G8prPTY+uBvir6/LeOYjxdKC4dPjhN/wkYK+lDFO+/pxu2JWmVpPdFMSZyp+MAbufSH8S9wL0ptE6heNPPj4hPNvGQkRRh9RmKN/3AiLhE0jSKsQYeUvHTQ/+vyhLeqph+m9a/Zyofv7Zifm3FugQcHU0PrN9WNgH+GRGDm7n/jSqXq3w+amaZFtcTEV9W8WXqSKBe0tCIeLVikZUU/wk1p8XXJCL+IWkQMAL4MsUPApzdXD0boPH2e7S0cET8KnXvjATukPSliGj4D7Qb1b8XOxz3AbdjKn4fq19F02CKboAFQG8VX9IhadPUz7YJRRfBPcB5FB8lt5C0S0TMjYjvU4wA9pFGm3qA1BcnaVeKo8/WhOFyip/V2VB/AL6qdGgm6WNNLHM/MCr1b/ah+ERQtSjG2n1O0rFpG0ohtUHLNTIdOE3pTBBJ27S0nvR6PBIRF1Ac6e7QaH1P8c6ng1ZT8WXtJhHxa2AcxZe4LT2vu4AzUnsXSVtRvCeOkNRTxQhpR6a2JkXEP4F/piNnqOjbVTG63V8jYjzFSGK7p/YPUHzHsXpDn+vGzgHcvm0BXK/0pRrpN8Oi+PmcY4Dvpy88ZlN8vOwC3JCOlP8CjE9/GGenL1LmUIzy9LtG2/kpsEl63GTg1PQRt1pXAb9X+hJuA/w3Rb/2nNRN8d9NLHMLxWhUjwMTKD7WttYJwBfSPptP8z8/VO1yAETE7yk+ms9U8eOrDV0pza3nhyq+HJtH0Z/7WKP1vQE8K6nqszca2Z7iE9NsirNTzl9PPaOB/dLrX0/xu3izKPq8H6X4pY6rI+Iv69nuacAVabuVnw6OA+al9o9SvH5Q/Cc6bcOeYsfg0dDM2iFJRwJDI2Jc7lrKkvqVx0TEU7lrycV9wGbtUETckj6id0gqBsCf2pnDF3wEbGaWjfuAzcwycQCbmWXiADYzy8QBbGaWiQPYzCyT/w/g8RbYojyH3QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "random.seed(1)\n",
    "perm_diffs = [perm_fun(session_times.Time, nA, nB) for _ in range(1000)]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.hist(perm_diffs, bins=11, rwidth=0.9)\n",
    "ax.axvline(x = mean_b - mean_a, color='black', lw=2)\n",
    "ax.text(50, 190, 'Observed\\ndifference', bbox={'facecolor':'white'})\n",
    "ax.set_xlabel('Session time differences (in seconds)')\n",
    "ax.set_ylabel('Frequency')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.121\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(perm_diffs > mean_b - mean_a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical Significance and P-Values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed difference: 0.0368%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeI0lEQVR4nO3de7xVZb3v8c83UBCVTEEjRZYZ2IYEgoXVxgtqioqJikoX8VJtOpqF5XGDQqW9ouPeXUzbplH6UowEwUQTdicBMxGV1iJUwGsJCpKuLaagHAH9nT/GWDjBBWsuZMxnXb7v12u+5hjPHGPM31iwvjw8c4xnKiIwM7PK+0DqAszM2ioHsJlZIg5gM7NEHMBmZok4gM3MEnEAm5klUlgAS+ou6T5JyyQtlTQmb79C0ipJi/PHSSX7XCbpWUlPSRpaVG1mZs2BiroOWFI3oFtELJK0J1ALnAqcBayLiB9vtX1v4DbgMOAjwBygV0S8XUiBZmaJtS/qwBGxGlidL6+V9ASw/3Z2GQ5MjYi3gOckPUsWxg9ta4cuXbpEVVXVzivarEC1tbUADBw4MHElVmm1tbX/ExFdt24vLIBLSaoCPgk8AgwGLpJ0DlADXBIRr5KF88Mlu61k+4FNVVUVNTU1hdRstrNJAvDf2TZI0oqG2gv/EE7SHsAdwMUR8TpwPXAw0J+sh/yTJh5vtKQaSTV1dXU7u1wzs4opNIAl7UIWvlMi4ncAEfFSRLwdEe8AvyIbZgBYBXQv2f2AvG0LETEpIqojorpr1/f06M3MWowir4IQcCPwRET8tKS9W8lmpwFL8uW7gc9L6iDpIKAnsLCo+szMUityDHgwMAp4XNLivO1y4AuS+gMBLAe+BhARSyXdDiwDNgFf9xUQZtaaFXkVxHxADbw0ezv7TAQmFlWTmVlz4jvhzMwScQCbmSXiADYzS8QBbGaWiAPYzCwRB7CZWSIOYDOzRBzAZmaJVGQ2NLOUqsbNKuzYy68aVtixrfVzD9jMLBEHsJlZIg5gM7NEHMBmZok4gM3MEnEAm5kl4gA2M0vEAWxmlogD2MwsEQewmVkiDmAzs0QcwGZmiTiAzcwScQCbmSXiADYzS8QBbGaWiAPYzCwRB7CZWSIOYDOzRBzAZmaJOIDNzBJxAJuZJeIANjNLxAFs1oxUVVUhqU08qqqqUv+4k2ufugAze9eKFSuIiNRlVISk1CUk5x6wmVkiDmAzs0QcwGYtwMqVKxk+fDg9e/bk4IMPZsyYMWzYsIGbb76Ziy66KHV577HHHnukLqFFcACbNXMRwemnn86pp57KM888w9NPP826desYP358Ie+3adOmQo5r7+UANmvm5s2bR8eOHTn//PMBaNeuHVdffTU33XQTb775Ji+88AJDhgyhZ8+eXHnllQC88cYbDBs2jH79+vGJT3yCadOmAVBbW8tRRx3FwIEDGTp0KKtXrwZgyJAhXHzxxVRXVzNx4kR69OjBO++8s/lY3bt3Z+PGjfztb3/jhBNOYODAgRxxxBE8+eSTADz33HN85jOf4dBDD2XChAmV/hG1WL4KwqyZW7p0KQMHDtyirXPnzhx44IFs2rSJhQsXsmTJEjp16sSgQYMYNmwYK1as4CMf+QizZs0C4LXXXmPjxo184xvf4K677qJr165MmzaN8ePHc9NNNwGwYcMGampqAFi0aBH3338/Rx99NPfccw9Dhw5ll112YfTo0dxwww307NmTRx55hAsvvJB58+YxZswYLrjgAs455xyuu+66yv6AWjAHsFkLd9xxx7HPPvsAcPrppzN//nxOOukkLrnkEsaOHcvJJ5/MEUccwZIlS1iyZAnHHXccAG+//TbdunXbfJyRI0dusTxt2jSOPvpopk6dyoUXXsi6detYsGABZ5555ubt3nrrLQAefPBB7rjjDgBGjRrF2LFjCz/v1sABbNbM9e7dmxkzZmzR9vrrr/P888/Tvn3791xPK4levXqxaNEiZs+ezYQJEzj22GM57bTT6NOnDw899FCD77P77rtvXj7llFO4/PLLWbNmDbW1tRxzzDG88cYb7LXXXixevLjB/X1db9N5DNismTv22GN58803mTx5MpD1XC+55BLOO+88OnXqxL333suaNWtYv349M2fOZPDgwbz44ot06tSJs88+m0svvZRFixZxyCGHUFdXtzmAN27cyNKlSxt8zz322INBgwYxZswYTj75ZNq1a0fnzp056KCDmD59OpB9OPjoo48CMHjwYKZOnQrAlClTiv6RtBoOYLNmThJ33nkn06dPp2fPnvTq1YuOHTvywx/+EIDDDjuMESNG0LdvX0aMGEF1dTWPP/44hx12GP379+fKK69kwoQJ7LrrrsyYMYOxY8fSr18/+vfvz4IFC7b5viNHjuQ3v/nNFkMTU6ZM4cYbb6Rfv3706dOHu+66C4BrrrmG6667jkMPPZRVq1YV+wNpRdSSb3usrq6O+g8NzLalatyswo69/KphZW9b/1/07f3OSWpTtyK3oXOtjYjqrdvdAzYzS8QBbGaWiAPYzCwRB7BZC3PFFVfw4x//mO9+97vMmTMHgAceeIA+ffrQv39/1q9fz6WXXkqfPn249NJLE1dr2+PrgM1aqO9///ubl6dMmcJll13G2WefDcCkSZNYs2YN7dq1K+tYmzZton17x0GluQds1gJMnDiRXr16cfjhh/PUU08BcN555zFjxgx+/etfc/vtt/Od73yHL33pS5xyyimsW7eOgQMHMm3aNOrq6hgxYgSDBg1i0KBBPPjgg0DWkx41ahSDBw9m1KhR293uy1/+MkOGDOGjH/0o11577ea6Jk+eTN++fenXrx+jRo0C2OZx7L38T55ZM1dbW8vUqVNZvHgxmzZtYsCAAVvMDfHVr36V+fPnc/LJJ3PGGWcA2Y0U9XesffGLX+Rb3/oWhx9+OM8//zxDhw7liSeeAGDZsmXMnz+f3XbbbbvbPfnkk9x3332sXbuWQw45hAsuuICnn36aH/zgByxYsIAuXbqwZs0aAMaMGbPN49iWCgtgSd2BycB+QACTIuIaSXsD04AqYDlwVkS8quwiyWuAk4A3gfMiYlFR9Zm1FA888ACnnXYanTp1ArLbhJtizpw5LFu2bPP666+/zrp16zYfa7fddmt0u2HDhtGhQwc6dOjAvvvuy0svvcS8efM488wz6dKlCwB77733do/jOYLfq8ge8CbgkohYJGlPoFbSvcB5wNyIuErSOGAcMBY4EeiZPz4FXJ8/m9n78M477/Dwww/TsWPH97xWOv/D9rbr0KHD5uV27dptd87g7R3HtlTYGHBErK7vwUbEWuAJYH9gOHBLvtktwKn58nBgcmQeBvaS1A2zNu7II49k5syZrF+/nrVr1/L73/++Sfsff/zx/PznP9+8vq3JdMrdrt4xxxzD9OnTeeWVVwA2D0E09ThtWUU+hJNUBXwSeATYLyJW5y/9g2yIArJwfqFkt5V529bHGi2pRlJNXV1dcUWbNRMDBgxg5MiR9OvXjxNPPJFBgwY1af9rr72Wmpoa+vbtS+/evbnhhhve13b1+vTpw/jx4znqqKPo168f3/72t3foOG1Z4XNBSNoDuB+YGBG/k/TPiNir5PVXI+JDku4BroqI+Xn7XGBsRGxzsgfPBWHl8FwQzVMbO9fKzwUhaRfgDmBKRPwub36pfmghf345b18FdC/Z/YC8zcysVSosgPOrGm4EnoiIn5a8dDdwbr58LnBXSfs5ynwaeK1kqMLMrNUp8iqIwcAo4HFJi/O2y4GrgNslfQVYAZyVvzab7BK0Z8kuQzu/wNrMzJIrLIDzsdxtfUfJsQ1sH8DXi6rHzKy58Z1wZs1Ijx492sx3q/Xo0SN1Cck5gM2akeXLl6cuwSrIk/GYmSXiADYzS8QBbGaWiAPYzCwRB7CZWSIOYDOzRBzAZmaJOIDNzBJxAJuZJeIANjNLxAFsZpaIA9jMLBEHsJlZIp4NzWwnaMr3zjVl26Z855y1PO4Bm5kl4gA2M0vEAWxmlogD2MwsEQewmVkiDmAzs0QcwGZmiTiAzcwScQCbmSXiADYzS8QBbGaWiAPYzCwRB7CZWSIOYDOzRBzAZmaJOIDNzBJxAJuZJeIANjNLxAFsZpaIA9jMLBEHsJlZIg5gM7NEHMBmZok4gM3MEnEAm5kl4gA2M0vEAWxmlogD2MwsEQewmVki7VMXYG1b1bhZhR17+VXDCju22c7gHrCZWSIOYDOzRBzAZmaJOIDNzBJxAJuZJeIANjNLpLAAlnSTpJclLSlpu0LSKkmL88dJJa9dJulZSU9JGlpUXWZmzUWRPeCbgRMaaL86Ivrnj9kAknoDnwf65Pv8QlK7AmszM0uusACOiD8Da8rcfDgwNSLeiojngGeBw4qqzcysOUgxBnyRpMfyIYoP5W37Ay+UbLMyb3sPSaMl1UiqqaurK7pWM7PCVDqArwcOBvoDq4GfNPUAETEpIqojorpr1647uTwzs8qpaABHxEsR8XZEvAP8ineHGVYB3Us2PSBvMzNrtSoawJK6layeBtRfIXE38HlJHSQdBPQEFlayNjOzSitsNjRJtwFDgC6SVgLfA4ZI6g8EsBz4GkBELJV0O7AM2AR8PSLeLqo2M7PmoLAAjogvNNB843a2nwhMLKoeM7PmpqwhCEmHFl2ImVlbU+4Y8C8kLZR0oaQPFlqRmVkbUVYAR8QRwJfIrlSolfRbSccVWpmZWStX9lUQEfEMMAEYCxwFXCvpSUmnF1WcmVlrVu4YcF9JVwNPAMcAn4uIf8mXry6wPjOzVqvcqyB+DvwauDwi1tc3RsSLkiYUUpmZWStXbgAPA9bXX5sr6QNAx4h4MyJuLaw6M7NWrNwx4DnAbiXrnfI2MzPbQeUGcMeIWFe/ki93KqYkM7O2odwAfkPSgPoVSQOB9dvZ3szMGlHuGPDFwHRJLwICPgyMLKooM7O2oKwAjoi/SPo4cEje9FREbCyuLDOz1q8pk/EMAqryfQZIIiImF1KVmVkbUFYAS7qV7JssFgP100QG4AA2M9tB5faAq4HeERFFFmNm1paUexXEErIP3szMbCcptwfcBVgmaSHwVn1jRJxSSFVmZm1AuQF8RZFFmJm1ReVehna/pB5Az4iYI6kT0K7Y0szMWrdyp6P8N2AG8Mu8aX9gZkE1mZm1CeV+CPd1YDDwOmyenH3foooyM2sLyg3gtyJiQ/2KpPZk1wGbmdkOKjeA75d0ObBb/l1w04HfF1eWmVnrV24AjwPqgMeBrwGzyb4fzszMdlC5V0G8A/wqf5hZAlXjZhVy3OVXDSvkuNa4cueCeI4Gxnwj4qM7vSIzszaiKXNB1OsInAnsvfPLMTNrO8oaA46IV0oeqyLiZ2Rf1GlmZjuo3CGIASWrHyDrETdlLmEzM9tKuSH6k5LlTcBy4KydXo2ZWRtS7lUQRxddiJlZW1PuEMS3t/d6RPx055RjZtZ2NOUqiEHA3fn654CFwDNFFGVm1haUG8AHAAMiYi2ApCuAWRFxdlGFmZm1duXeirwfsKFkfUPeZmZmO6jcHvBkYKGkO/P1U4FbCqnIzKyNKPcqiImS/hs4Im86PyL+WlxZZmatX7lDEACdgNcj4hpgpaSDCqrJzKxNKPcrib4HjAUuy5t2AX5TVFFmZm1BuT3g04BTgDcAIuJFYM+iijIzawvKDeANERHkU1JK2r24kszM2oZyA/h2Sb8E9sq/IXkOnpzdzOx9afQqCEkCpgEfJ/tW5EOA70bEvQXXZmbWqjUawBERkmZHxKGAQ9fMbCcpdwhikaRBhVZiZtbGlHsn3KeAsyUtJ7sSQmSd475FFWZm1tptN4AlHRgRzwNDK1SPmVmb0VgPeCbZLGgrJN0RESMqUJOZWZvQ2BiwSpb9FfRmZjtRYwEc21g2M7P3qbEhiH6SXifrCe+WL8O7H8J1LrQ6M7NWbLs94IhoFxGdI2LPiGifL9evbzd8Jd0k6WVJS0ra9pZ0r6Rn8ucP5e2SdK2kZyU9JmnAzjk9M7PmqynTUTbVzcAJW7WNA+ZGRE9gbr4OcCLQM3+MBq4vsC4zs2ah3OuAmywi/iypaqvm4cCQfPkW4E9k01wOBybnE/48LGkvSd0iYnVR9VnjqsbNKuzYy68aVtixzVqKInvADdmvJFT/wbvfK7c/8ELJdivztveQNFpSjaSaurq64io1MytYpQN4s9LpLZu436SIqI6I6q5duxZQmZlZZVQ6gF+S1A0gf345b18FdC/Z7oC8zcys1ap0AN8NnJsvnwvcVdJ+Tn41xKeB1zz+a2atXWEfwkm6jewDty6SVgLfA64im9z9K8AK4Kx889nAScCzwJvA+UXVZWbWXBR5FcQXtvHSsQ1sG8DXi6rFzKw5SvYhnJlZW+cANjNLxAFsZpaIA9jMLBEHsJlZIg5gM7NEHMBmZok4gM3MEnEAm5kl4gA2M0vEAWxmlogD2MwsEQewmVkiDmAzs0QcwGZmiTiAzcwScQCbmSXiADYzS8QBbGaWiAPYzCwRB7CZWSIOYDOzRBzAZmaJOIDNzBJxAJuZJeIANjNLxAFsZpaIA9jMLBEHsJlZIg5gM7NEHMBmZok4gM3MEnEAm5kl4gA2M0vEAWxmlogD2MwsEQewmVki7VMXYGbNR9W4WYUcd/lVwwo5bkvnHrCZWSLuAbdg7q2YtWzuAZuZJeIANjNLxAFsZpaIA9jMLBEHsJlZIg5gM7NEHMBmZok4gM3MEnEAm5kl4gA2M0vEAWxmlogD2MwskSST8UhaDqwF3gY2RUS1pL2BaUAVsBw4KyJeTVGfmVklpOwBHx0R/SOiOl8fB8yNiJ7A3HzdzKzVak5DEMOBW/LlW4BT05ViZla8VAEcwB8l1UoanbftFxGr8+V/APs1tKOk0ZJqJNXU1dVVolYzs0KkmpD98IhYJWlf4F5JT5a+GBEhKRraMSImAZMAqqurG9zGzKwlSNIDjohV+fPLwJ3AYcBLkroB5M8vp6jNzKxSKh7AknaXtGf9MnA8sAS4Gzg33+xc4K5K12ZmVkkphiD2A+6UVP/+v42IP0j6C3C7pK8AK4CzEtRmZlYxFQ/giPg70K+B9leAYytdj5lZKs3pMjQzszbFAWxmlogD2MwsEQewmVkiDmAzs0QcwGZmiTiAzcwScQCbmSXiADYzS8QBbGaWiAPYzCwRB7CZWSIOYDOzRBzAZmaJOIDNzBJxAJuZJeIANjNLxAFsZpZIqq+lN7M2qmrcrEKOu/yqYYUct0juAZuZJeIANjNLxAFsZpaIA9jMLBEHsJlZIg5gM7NEHMBmZok4gM3MEnEAm5kl4gA2M0vEtyIXxLdbmllj3AM2M0vEAWxmlogD2MwsEQewmVkiDmAzs0QcwGZmiTiAzcwScQCbmSXiADYzS8QBbGaWiAPYzCwRB7CZWSIOYDOzRBzAZmaJOIDNzBLxfMBm1uq0lPm422QAt5Q/HDNr3TwEYWaWiAPYzCwRB7CZWSIOYDOzRJpdAEs6QdJTkp6VNC51PWZmRWlWASypHXAdcCLQG/iCpN5pqzIzK0azCmDgMODZiPh7RGwApgLDE9dkZlaI5hbA+wMvlKyvzNvMzFodRUTqGjaTdAZwQkR8NV8fBXwqIi4q2WY0MDpfPQR4quKFvj9dgP9JXUQBfF4tS2s8r+Z8Tj0iouvWjc3tTrhVQPeS9QPyts0iYhIwqZJF7UySaiKiOnUdO5vPq2VpjefVEs+puQ1B/AXoKekgSbsCnwfuTlyTmVkhmlUPOCI2SboI+L9AO+CmiFiauCwzs0I0qwAGiIjZwOzUdRSoxQ6fNMLn1bK0xvNqcefUrD6EMzNrS5rbGLCZWZvhAC6YpL0l3Svpmfz5Qw1s01/SQ5KWSnpM0sgUtTZFOeeVb/cHSf+UdE+la2yKxm6Bl9RB0rT89UckVSUos0nKOKcjJS2StCm/BLRFKOO8vi1pWf67NFdSjxR1lsMBXLxxwNyI6AnMzde39iZwTkT0AU4AfiZpr8qVuEPKOS+AHwGjKlbVDijzFvivAK9GxMeAq4H/qGyVTVPmOT0PnAf8trLV7bgyz+uvQHVE9AVmAP9Z2SrL5wAu3nDglnz5FuDUrTeIiKcj4pl8+UXgZeA9F203M42eF0BEzAXWVqimHVXOLfCl5zsDOFaSKlhjUzV6ThGxPCIeA95JUeAOKue87ouIN/PVh8nuJ2iWHMDF2y8iVufL/wD2297Gkg4DdgX+VnRh71OTzquZK+cW+M3bRMQm4DVgn4pUt2Na6239TT2vrwD/XWhF70OzuwytJZI0B/hwAy+NL12JiJC0zctOJHUDbgXOjYjkvZKddV5mKUg6G6gGjkpdy7Y4gHeCiPjstl6T9JKkbhGxOg/Yl7exXWdgFjA+Ih4uqNQm2Rnn1UI0egt8yTYrJbUHPgi8Upnydkg559QSlXVekj5L1lE4KiLeqlBtTeYhiOLdDZybL58L3LX1Bvlt13cCkyNiRgVrez8aPa8WpJxb4EvP9wxgXjTvi+hb6239jZ6XpE8CvwROiYjm3TGICD8KfJCNE84FngHmAHvn7dXAr/Pls4GNwOKSR//Utb/f88rXHwDqgPVk43VDU9e+jfM5CXiabOx9fN72fbJfYoCOwHTgWWAh8NHUNe+EcxqU/5m8QdabX5q65p10XnOAl0p+l+5OXfO2Hr4TzswsEQ9BmJkl4gA2M0vEAWxmlogD2MwsEQewmVkiDmB7D0kfljRV0t8k1UqaLalX6rrqSTqloVmwCn7PiyV12oH9fibpyCJq2lGSTi2dwEbSjyUdk7KmtsqXodkW8glmFgC3RMQNeVs/oHNEPFDhWtpFxNsVei+R/T40eAu4pOVkM2yV/a27kvYBZkXEp3dCfe0jm4PifZN0M3BP5Df95NM1/ioijt8Zx7fyuQdsWzsa2FgfvgAR8WhEPKDMjyQtkfR4/bzFkoZI+pOkGZKelDQl3/YESdPrj5Nvd0++fHw+B/IiSdMl7ZG3L5f0H5IWAWdK+mbJ3K5T823Ok/Rf+XKVpHklc78emLffLOlaSQsk/b2h+W7zfZ+SNBlYAnSXdL2kGmVzM1+Zb/dN4CPAfZLu2179WxkB/KHk/ZZL+s/8Z7dQ0sfy9q6S7pD0l/wxOG+/QtKtkh4EbpW0n6Q7JT2aP/413+7s/HiLJf1S2ZSNSFonaWK+7cP5/v8KnAL8KN/+4IhYAewjqaF5P6xIqe8E8aN5PYBvAldv47URwL1kX5i6H9l8st2AIWSzgx1A9o/6Q8DhZHONPA/snu9/Pdldf12AP5e0jwW+my8vB/695D1fBDrky3vlz+cB/5Uv/55s8iKALwMz8+Wbye5c+wDZvLHPNnA+VWRTMX66pK3+jr52wJ+AviV1dcmXt1n/Vse/Bfhcyfpy3r1z6xyyXihk8/Eeni8fCDyRL18B1AK75evTgItL6vsg8C/5z2CXvP0XZHNLA0T9+5PNiTuh5Gdzxla1/goYkfrvX1t7eDIea4rDgdsiGxZ4SdL9ZLezvg4sjIiVAJIWA1URMV/SH4DPSZoBDAP+nWx2qt7Ag9n//NmVLLTrTStZfgyYImkmMLOBmj4DnJ4v38qWk2/PjGxIYZmkbU2XuSK2nPzoLEmjyf7x6JbX+dhW+3y6kfrrdSO7DbvUbSXPV+fLnwV6693phTuX9Kjvjoj1+fIxZMFN/mfwmqRRwEDgL/n+u/HuxEgbgPpvIqkFjmugxnovk/XyrYIcwLa1pWSTzTRV6YxTb/Pu362pwEXAGqAmItbm4633RsQXtnGsN0qWhwFHAp8Dxks6dAdr2tbk6ZvfS9JBwP8GBkXEq/lYaccG9mms/nrrG9g/Glj+AFkv/P9t8SZZoJb+LBoisvH6yxp4bWPk3Vu2/DNpSMe8XqsgjwHb1uYBHfJeIACS+ko6gmxinZGS2knqShaMCxs53v3AAODfyMIYsm8pGFwyBrq7GrjKQtIHgO4RcR/Zf/M/CGw91rqAbEYsgC/lNe6ozmSB91reYz6x5LW1wJ5NqR94AvjYVm0jS57re81/BL5Rv4Gk/tuoby5wQb5NO0kfzNvOkLRv3r63Gv8OtNJzqdeLbBzcKsgBbFvIe0ynAZ9VdhnaUuD/kH3rxZ1k/x1/lCyo/z0i/tHI8d4m+2/wifkzEVFHNo57m6THyILo4w3s3g74jaTHyb7n69qI+OdW23wDOD8/zihgTFPPuaTWR/P3eZJsXPbBkpcnAX+QdF8T6p9FNj5e6kP5PmOAb+Vt3wSq8w8SlwH/axsljgGOzn8etUDviFgGTAD+mB/3XrKhj+2ZClwq6a+SDpa0C9k/FDWN7Gc7mS9DMyuQpPnAyRHxT+3ApWyVIOk0YEBEfCd1LW2Ne8BmxbqE7MqG5qw98JPURbRF7gGbmSXiHrCZWSIOYDOzRBzAZmaJOIDNzBJxAJuZJeIANjNL5P8DGAhl59HtJIQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "random.seed(1)\n",
    "obs_pct_diff = 100 * (200 / 23739 - 182 / 22588)\n",
    "print(f'Observed difference: {obs_pct_diff:.4f}%')\n",
    "conversion = [0] * 45945\n",
    "conversion.extend([1] * 382)\n",
    "conversion = pd.Series(conversion)\n",
    "\n",
    "perm_diffs = [100 * perm_fun(conversion, 23739, 22588) \n",
    "              for _ in range(1000)]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.hist(perm_diffs, bins=11, rwidth=0.9)\n",
    "ax.axvline(x=obs_pct_diff, color='black', lw=2)\n",
    "ax.text(0.06, 200, 'Observed\\ndifference', bbox={'facecolor':'white'})\n",
    "ax.set_xlabel('Conversion rate (percent)')\n",
    "ax.set_ylabel('Frequency')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## P-Value\n",
    "If `np.mean` is applied to a list of booleans, it gives the percentage of how often True was found in the list (#True / #Total)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.332\n"
     ]
    }
   ],
   "source": [
    "print(np.mean([diff > obs_pct_diff for diff in perm_diffs]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-value for single sided test: 0.3498\n"
     ]
    }
   ],
   "source": [
    "survivors = np.array([[200, 23739 - 200], [182, 22588 - 182]])\n",
    "chi2, p_value, df, _ = stats.chi2_contingency(survivors)\n",
    "\n",
    "print(f'p-value for single sided test: {p_value / 2:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# t-Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-value for single sided test: 0.1408\n"
     ]
    }
   ],
   "source": [
    "res = stats.ttest_ind(session_times[session_times.Page == 'Page A'].Time, \n",
    "                      session_times[session_times.Page == 'Page B'].Time,\n",
    "                      equal_var=False)\n",
    "print(f'p-value for single sided test: {res.pvalue / 2:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-value: 0.1408\n"
     ]
    }
   ],
   "source": [
    "tstat, pvalue, df = sm.stats.ttest_ind(\n",
    "    session_times[session_times.Page == 'Page A'].Time, \n",
    "    session_times[session_times.Page == 'Page B'].Time,\n",
    "    usevar='unequal', alternative='smaller')\n",
    "print(f'p-value: {pvalue:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ANOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAENCAYAAADUlXqkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZM0lEQVR4nO3df5RdZX3v8feHJCIFwUqAokgmFrBDBsgllOJtwIzpxdZf9LZWmXVRrk7rpRfS5dVyiWuo/Fhr7sVeL3eJqEgdBLt0EI26EBRFc0YIq0iTS8DAyDX8kiCVApZfppAM3/vH3lNOJmdm9pmcZ59z9nxea+2Vc56zf3yfZPKdfZ699/NVRGBmlsJe7Q7AzKrLCcbMknGCMbNknGDMLBknGDNLZmG7A9gTixcvjp6enlKP+fzzz7PvvvuWeswyVb1/UP0+tqN/mzZteiIiDpra3tUJpqenh40bN5Z6zLGxMVatWlXqMctU9f5B9fvYjv5JerhRu78imVkyTjBmlowTjJkl4wRjZsk4wZhZMskSjKSrJD0uaUtd23GS/kHSTyR9W9L+eXuPpO2SNufLFaniMrPypDyDuRr4wyltXwDWRsQxwDeBc+s+uz8ilufLWQnjMrOSJEswEXEL8NSU5qOAW/LXNwN/mur4ZtZ+Zd9odw9wGvAt4M+A19d9tlTSncAzwPkRcWujHUj6EPAhgEMOOYSxsbGWB9nf3z+n7Wq1WosjKd9zzz2X5O+0k1S9jx3Vv4hItgA9wJa6978DfB/YBFwAPJm37w0cmL9eATwC7D/b/lesWBFlW3LeDaUfs0y1Wq3dISRX9T62o3/Axmjwf7TUM5iI+ClwKoCko4C35+0vAC/krzdJup/s61S5zwGYWUuVepla0sH5n3sB5wNX5O8PkrQgf/0G4EjggTJjM7PWS3YGI2kUWAUslrSN7CvRfpLOzlf5BvDF/PUpwMWSdgAvAWdFxNQBYjPrMskSTEQMTPPRpxqsuw5YlyoWM2sP38lrZsk4wZhZMk4wNm+Mjo7S19fH6tWr6evrY3R0tN0hVV5Xz2hnVtTo6ChDQ0OMjIwwMTHBggULGBwcBGBgYLrhQttTPoOxeWF4eJiRkRH6+/tZuHAh/f39jIyMMDw83O7QKs0JxuaF8fFxVq5cuUvbypUrGR8fb1NE84MTjM0Lvb29bNiwYZe2DRs20Nvb26aI5gcnGJsXhoaGGBwcpFarsXPnTmq1GoODgwwNDbU7tErzIK/NC5MDuWvWrGF8fJze3l6Gh4c9wJuYE4zNGwMDAwwMDFS+LlIn8VckM0vGCcbMknGCMbNknGDMLBknGDNLxleR5ilJc9oum37VrBifwcxTjSZojgiWnHfDbBO5mxXWEZUd888+JmmrpPskvTVVXGZWno6o7CjpaOB0YFm+zWcnJwE3s+7VKZUdTwOujYgXIuJBYCtwYqrYzKwcnVLZ8XXA7XXrbcvbdlNGZcfZdEzVvESq3r+OqnyYQCf1r+wE80HgMkl/A1wPvNjsDiLiSuBKgBNOOCFKf6bkphur/RxL1fsHlX8WqZP61xGVHYFH2bVO9WF5m5l1sY6o7Eh2NnO6pL0lLSWr7HhHmbGZWet1RGXHiLhH0nXAvcBO4OyImEgVm5mVoyMqO+brDwOegdmsQnwnr5kl4wRjZsk4wZhZMk4wZpaME4yZJeMEY2bJOMGYWTJOMGaWjBOMmSXjBGNmyTjBmFkyTjBmlowTjJkl4wRjZsk4wZhZMk4wZpZM2YXXlku6XdJmSRslnZi3r5L0dN6+WdLHU8VlZuUpu/Da3wIXRcRy4OP5+0m3RsTyfLk4YVxmVpKyC68FMFku9gDgF6mOb2btV3ZdpA8D35P0SbLk9u/rPnuTpLvIks5fR8Q9jXbgwmvFnf3D53l+R/Pb9ay9san1910En1m9b/MHapNOKkyWQkf1LyKSLUAPsKXu/WXAn+av3wP8IH+9P7Bf/vptwM+K7H/FihVRtiXn3VD6MedqLrHWarVSjtNOc+ljN2lH/4CN0eD/aNlXkc4kK1cC8DXy+tMR8UxEPJe//g6wSNLikmMzsxYrO8H8Anhz/votwM8AJP2WJOWvT8zjerLk2MysxcouvPYXwKckLQT+lXwsBXg38JeSdgLbgdPz0y4z62LtKLy2osG6lwOXp4rFzNrDd/KaWTJOMGaWjBOMmSXjBGNmycw6yCvplcA7gJOB15Jd5dkC3BjT3G1rZgazJBhJF5EllzHgx8DjwCuBo4BL8uTz0Yi4O3GcZtaFZjuDuSMiLpjms0slHQwc3uKYzKwiZkwwEbHbU2+S9iJ7buiZiHic7KzGzGw3hQZ5JX1F0v6S9iUbf7lX0rlpQzOzblf0KtLREfEM8MfAd4GlwPtSBWVm1VD0UYFFkhaRJZjLI2KHJD8r1OFe1buWY65Z2/yG1zR7HIC3N38cq7yiCebzwEPAXcAtkpYAz6QKylrj2fFLeOiS5v7jj42NsWrVqqa2aXaCKps/CiWYiLiMbLKoSQ9L6k8TkplVxWz3wXxklu0vbWEsZlYxs53BvCr/843A7wLX5+/fCdyRKigzq4bZ7oO5CEDSLcDxEfFs/v5CwF+8zWxGRS9THwK8WPf+xbzNzGxaRRPMl4A7JF2Yn738mAIXM5us7ihJl0naKuluScfPoT9m1kEKJZiIGAY+CPwqXz4QEf+jwKZXU7y64x8BR+bLh4DPFYnNzDpXM3PybgYem9xG0uER8fOZNoiIWyT1TG2mcXXH04Av5ZN93y7p1ZIOjYjHmojRzDpIoQQjaQ1ZVYBfAhOAyBLFsXM45odpXN3xdcAjdetty9t2STCu7NicZmOda1XAbvo76ajKhwl0VP8aVWObugBbgQOLrNtg2x6KVXe8AVhZt94PgRNm2rcrO87MlR0bc2XH1mMPKzs+AjzdopzWsLoj8Cjw+rr1DsvbzKxLFR2DeQAYk3Qj8MJkY0TM5U7eyeqOY9RVdyS7ie8cSdcCvwc8HR5/MetqRRPMz/PlFflSSJPVHb9DVvh+K/Br4ANFj2Nmnanow46Td/Tul79/ruB2zVR3DODsIvs1s+5QdEa7Pkl3AvcA90jaJGlZ2tDMrNsVHeS9EvhIRCyJiCXAR4G/SxeWmVVB0QSzb0TUJt9ExBiwb5KIzKwyCl9FkvQ3wN/n788gu7JkZjatognmg8BFZPevBHBr3mYdbk7TWd7U3DYH7LOo+WPYvFD0KtKvgL9KHIu1WLPz8UKWkOaynVkjRa8i3Szp1XXvf1PS95JFZWaVUHSQd3FE/Mvkm/yM5uAkEZlZZRRNMC9J+rca1HnZEtdFMrMZFR3kHQI2SPoR2VQNJ/PyLf5d67iLvs/T23c0vV0zA6cH7LOIuy44teljpCZp+s8+Mf122Q3XZsUUHeS9KZ/C8qS86cMR8US6sMrx9PYdyQuTdWpRsukSxVwKr5lNp+ggr8imvjw+Im4AfmNyLl0zs+kUHYP5LPAmYPLhxWeBzySJyMwqo+gYzO9FxPH5A49ExK8kFZ62wczmp6IJZoekBeRXjiQdBLyULCqzgqYbqH/4E++Y0/6WnHfDbm2dOlDfDYommMuAbwIHSxoG3g2cnywqs4KmHai/ZPqrXVUZqO8GRa8ifVnSJmA12WXqP46I8Zm2kXQV8A7g8Yjoy9u+SlbnGuDVwL9ExPK8tMk4cF/+2e0RcVaTfTGzDlO0bMlvAw9GxGckrQL+g6TH6u/ubeBq4HKyqpAARMR76/b5v9l1IvH7IyvGZmYVUfQq0jpgQtIRwOfJZv//ykwbRMQtwFONPssve78HGC0eqpl1m8KPCkTETuBPgMsj4lzg0D047snALyPiZ3VtSyXdKelHkk7eg32bWYdo5irSAPB+4J15255MAjLArmcvjwGHR8STklYA35K0LCKembphqys7llH5sGOq7BXQUVUBC3hV71qOuWZt8xte08wxYGyseyZw7Kh/w0bV2KYuwNFkV5IG8vdLgfMKbNdDXVXHvG0hWQnaw2bYboxZqjpGCyo7llH50FUP0/K/4e46qbJj0atI91I34VREPAjM8EjcjP4A+GlEbJtsyO+reSoiJiS9ATgST8lp1vWKjsE0LS+69g/AGyVtkzSYf3Q6uw/ungLcLWkz8HXgrIhoOEBsZt2j6BhM02KaomsR8Z8btK0ju1JlZhWS7AzGzKzojXZHAecCS+q3iYi3JIrLzCqg6FekrwFXkFVznEgXjplVSdEEszMiPpc0EjOrnKJjMN+W9F8lHSrpNZNL0sjMrOsVPYM5M//z3Lq2AN7Q2nDMrEqK3mi3NHUgZlY9MyYYSW+JiPWS/qTR5xHxjTRhmVkVzHYG82ZgPS8/4FgvACcYM5vWjAkmIi7I//xAOeGYWZXM9hXpDOArEdFwgu98prtDI2JDiuBSK+tRf2iuuJtZVcz2FelA4M58Pt5NwD8DrwSOIPv69AQwh/+hneHZ8UvmbWVHszLM9hXpU5IuB94C/D5wLLCdbILu90XEz9OHaGbdatbL1BExAdycL2ZmhflpajNLJtl8MGZlaTTO1erKjjY3TjDW1aYdpG9hZUebu0JfkSQdImlE0nfz90fXTYE53TZXSXpc0pa6tq9K2pwvD+VTZE5+9jFJWyXdJ+mtc+yPmXWQomMwVwPfA16bv/9/wIcLbPOH9Q0R8d6IWB5ZBcd15HcCSzqabK7eZfk2n5W0oGBsZtahiiaYxRFxHfASQGRF2GaceCqaq+x4GnBtRLyQVyzYCpxYMDYz61BFx2Cel3Qg2fNHSDqJXetKN2tqZcfXAbfXfb4tb9uNC6+l1VFFuxKpSh/7+/vntF2tVmtxJDNoVCxp6gIcD9xGllRuI/uKdGyB7XqYUngtb/8c8NG695cDZ9S9HwHePdv+XXit9bqt8NpcVL2P7fiZYw8Lr/1fSW8G3ggIuC8idswloUlaSFbjekVd86PA6+veH5a3mVkXK1pVYAHwNrIzkoXAqZKIiEvncMzdKjsC1wNfkXQp2UDykcAdc9i3mXWQomMw3wb+FfgJ+UDvbPLKjquAxZK2ARdExAgNKjtGxD2SrgPuBXYCZ0f2iIKZdbGiCeawiDi2mR1HE5Ud8/ZhYLiZY5hZZyt6mfq7kk5NGomZVU7RM5jbgW9K2gvYQTbQGxGxf7LIzKzrFU0wlwJvAn6SX5IyM5tV0a9Ij5Ddz+LkYmaFFT2DeQAYyx92fGGycY6Xqc2sCcdd9H2e3t7cbWfNTtV6wD6LuOuC1g+zFk0wD+bLK/LFzEry9PYdTc0dPZfpKFLNHV30Tt6LkhzdzCpttrIll0fEOZK+Tf6gY72IeFeyyMys6812BvN+4BzgkyXEYmYVM1uCuR8gIn5UQixmVjGzJZiDJH1kug99FcnMZjJbglkA7Ed2566ZWVNmSzCPRcTFpURiZpUz2528PnMxszmb7QxmdSlRmNm0XtW7lmOuWdvcRtc0ewyA4jfzFTVjgomIhlUBzKw8z45f0rV38ro2tZklkyzBNKrsmLevkfRTSfdI+tu8rUfS9rqqj1ekisvMypOyNvXVZOVIvjTZIKmfrMjacRHxgqSD69a/P7KKj2ZWEcnOYKJxZce/BC6JiBfydR5PdXwza7+UZzCNHAWcLGmYrErBX0fEP+afLZV0J/AMcH5E3NpoB67smFZVqh7OpBv72Ey8c+1fkr+TRtXYWrUwpbIjsAX4NNn9NSeSzTEjYG/gwHydFWQz6O0/2/5d2bH1ql71MKL7+tjsz9Bc+renP6dMU9mx7KtI24Bv5DHdQVZjaXFkRe+fBIiITWQPWR5Vcmxm1mJlJ5hvAf0Ako4imx3vCUkH5dUjkfQGssqOD5Qcm5m1WLIxmEaVHYGrgKvyS9cvAmdGREg6BbhY0g6ys5qzwjf5mXW9ZAkmpqnsCJzRYN11wLpUscxkTncw3lR8mwP2WdT8/s2maPRz+vAn3jGnfS0574bd2lL9nJZ9FamjNHP79aSetTfOaTuzuZr25+2SxlWE5vKoQCp+VMDmjdHRUfr6+li9ejV9fX2Mjo62O6TKm9dnMDZ/jI6OMjQ0xMjICBMTEyxYsIDBwUEABgam+zZve8pnMDYvDA8PMzIyQn9/PwsXLqS/v5+RkRGGh4fbHVqlOcHYvDA+Ps7KlSt3aVu5ciXj4+Ntimh+cIKxeaG3t5cNGzbs0rZhwwZ6e3vbFNH84ARj88LQ0BCDg4PUajV27txJrVZjcHCQoaGhdodWaR7ktXlhciB3zZo1jI+P09vby/DwsAd4E3OCsXljYGCAgYGBjrpPpOr8FcnMknGCMbNknGDMLBknGDNLxgnGzJJxgjGzZJxgzCyZjii8lrd/TNJWSfdJemuquMysPB1ReE3S0cDpwDLgtcAPJB0VERMJ4zOzxDql8NppwLV5dYEHga1kZU3MrIt1SuG11wG31623LW/bTasLr81FtxXtakY3FiVrVtX72En9KzvBLAReA5wE/C5wXV6mpLCIuBK4EuCEE06I0p8puenGSj/HMh+e06l6Hzupfx1ReA14FHh93XqH5W1m1sU6ovAacD1wuqS9JS0lK7x2R8mxmVmLdUThNeAeSdcB9wI7gbN9Bcms+3VE4bV8/WHAMzCbVYjv5DWzZJxgzCwZJxgzS8YJxsyScYIxs2ScYMwsGScYM0vGCcbMknGCMbNknGDMLBknGDNLxgnGzJJxgjGzZJxgzCyZsqfM7AqSZv78E43bs6ltzGySz2AaiIhpl1qtNu1nZrYrJxgzS6bUyo6SLpT0qKTN+fK2vL1H0va69itSxWVWVaOjo/T19bF69Wr6+voYHR1td0jlVnbM/Z+I+GSD9e+PiOUJ4zGrrNHRUYaGhhgZGWFiYoIFCxYwODgIwMDAdLPXpld2ZUczS2B4eJiRkRH6+/tZuHAh/f39jIyMMDzc3mmu23EV6RxJ7wc2Ah+NiF/l7Usl3Qk8A5wfEbc22rjdlR07qWpeClXvH1Szj+Pj40xMTDA2NvZv/ZuYmGB8fLy9fZ3pismeLkAPsKXu/SHAArIzp2Hgqrx9b+DA/PUK4BFg/9n2v2LFiihbrVYr/Zhlqnr/IqrZx2XLlsX69esj4uX+rV+/PpYtW1bK8YGN0eD/aKlXkSLilxExEREvAX9HXuA+sqL3T+avNwH3k9WxNrMChoaGGBwcpFarsXPnTmq1GoODgwwNDbU1rlK/Ikk6NCIey9/+R2BL3n4Q8FRETOS1qo8EHigzNrNuNjmQu2bNGsbHx+nt7WV4eLitA7xQfmXHVZKWAwE8BPyXfPVTgIsl7SCrV31WRHiA2KwJAwMDDAwMMDY2xqpVq9odDlB+ZceRadZdB6xLFYuZtYfv5DWzZJxgzCwZJxgzS8YJxsySUXTxNAOS/hl4uOTDLgaeKPmYZap6/6D6fWxH/5ZExEFTG7s6wbSDpI0RcUK740il6v2D6vexk/rnr0hmlowTjJkl4wTTvCvbHUBiVe8fVL+PHdM/j8GYWTI+gzGzZJxgzCyZeZdgJE3kE4tvkfQ1Sb+R8FjnSNoqKSQtTnWcKccss39flnRffqyrJC1Kdawpxy2zjyOS7pJ0t6SvS9ov1bGmHLe0PtYd8zJJz7Vyn/MuwQDbI2J5RPQBLwJnJTzWbcAfUO7NgGX278vA7wDHAPsAf57wWPXK7ON/i4jjIuJY4OfAOQmPVa/MPiLpBOA3W73f+Zhg6t0KHCHpnZJ+LOlOST+QdAhkE2FJulnSPZK+IOnhyTMRSWdIuiP/LfN5SQum7jwi7oyIh8rt0i5S9+87dVMm3gEcVmrvMqn7+Ey+rsiSaDuuiiTtY972v4D/3urA522CkbQQ+CPgJ8AG4KSI+HfAtbz8F30BsD4ilgFfBw7Pt+0F3gv8fmSlViaA/1RqB2ZRZv/yr0bvA25K0pnpj1tKHyV9EfgnsrO1T6fqzzTHLqOP5wDX18022TLzsTb1PpI2569vJZsE643AVyUdCrwCeDD/fCXZ1J5ExE2SJisgrCabnPwfs19s7AM8Xkr0s2tH/z4L3BLTVIJIoNQ+RsQH8t/ynyb7D/vFVneogVL6KOm1wJ+RzT7Zeo1mAq/yAjzXoG0MeFf+ehUwlr/eDCytW+8psgfJ1gD/s4ljPgQsrmL/yH57fgvYq8r/hvm2pwA3VKmPwNvJzs4eypeXgK2t6se8/Yo0xQHAo/nrM+vabwPeAyDpVF4eBPsh8G5JB+efvUbSkpJinYsk/ZP058BbgYHIKkW0U8v7qMwRk6+BdwE/TdaD2bW8jxFxY0T8VkT0REQP8OuIOKJlEZf1W6dTFhr/ZjiNrIrBJrLBrsnfDAfn/0hbyMqsPAbsnX/2XrLfHHfn253UYL9/BWwDdgK/AL5Qsf7tJCsxszlfPl6lf0OyMcrbyMY/tpBdNZu1Xlc39bHIcfdk8aMCM5C0NzARETslvQn4XFSofnbV+wfuY7vNx0HeZhwOXCdpL7J7Ef6izfG0WtX7B+5jW/kMxsyS8SCvmSXjBGNmyTjBmFkyHuS1lpM0QXZpdyEwDpwZEb9ub1TWDj6DsRRKfRLYOpcTjKWW9Elg62xOMJZM1Z9Yt9l5DMZSqPoT61aQE4ylsH3qreqSPg1cGhHXS1oFXDjLPgRcExEfSxGglcNfkawsVX9i3RpwgrGyXAh8TdImdi3MfhFwqqQtZBMf/RPwbETcC5wPfF/S3cDNwKHlhmx7ys8iWVt18pPAtuc8BmPt1rFPAtue8xmMmSXjMRgzS8YJxsyScYIxs2ScYMwsGScYM0vm/wN/ZTXPNIDmgAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "four_sessions = pd.read_csv(FOUR_SESSIONS_CSV)\n",
    "\n",
    "ax = four_sessions.boxplot(by='Page', column='Time',\n",
    "                           figsize=(4, 4))\n",
    "ax.set_xlabel('Page')\n",
    "ax.set_ylabel('Time (in seconds)')\n",
    "plt.suptitle('')\n",
    "plt.title('')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     Page  Time\n",
      "0  Page 1   164\n",
      "1  Page 2   178\n",
      "2  Page 3   175\n",
      "3  Page 4   155\n",
      "4  Page 1   172\n"
     ]
    }
   ],
   "source": [
    "print(pd.read_csv(FOUR_SESSIONS_CSV).head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed means: [172.8 182.6 175.6 164.6]\n",
      "Variance: 55.426666666666655\n",
      "30.146666666666615\n"
     ]
    }
   ],
   "source": [
    "observed_variance = four_sessions.groupby('Page').mean().var()[0]\n",
    "print('Observed means:', four_sessions.groupby('Page').mean().values.ravel())\n",
    "print('Variance:', observed_variance)\n",
    "# Permutation test example with stickiness\n",
    "def perm_test(df):\n",
    "    df = df.copy()\n",
    "    df['Time'] = np.random.permutation(df['Time'].values)\n",
    "    return df.groupby('Page').mean().var()[0]\n",
    "    \n",
    "print(perm_test(four_sessions))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pr(Prob) 0.087\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfDElEQVR4nO3dfZRU1Znv8e9jgyFIFMSWEECaSVqMCI10Q/A6RrSHiMAIKoi5RojLO7iiSTDhJhBlzehdJiFrzY2BuQ4JUUZQImAbhKskuYhoVnxNNwEBIQEjBJCXVrQNYAT0uX+c3Vh0gC7aPr2rqn+ftWrVOfu81FNV9s/DrnP2MXdHRERa3imxCxARaa0UwCIikSiARUQiUQCLiESiABYRiaRN7AI+jrPOOstLSkpilyF5rqamBoDy8vLIlUihqqmpedPdixu2Wz6fhlZRUeHV1dWxy5A8Z2YA5PPfguQ2M6tx94qG7eqCEBGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUjyejjKXFYy9clU9rtl+ohU9isiLU9HwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJJLUANrPeZrY64/Gumd1uZmea2XIz2xSeO4X1zcxmmtlmM3vFzAakVZuISC5ILYDd/Y/u3t/d+wPlwAFgMTAVWOHupcCKMA9wJVAaHhOBWWnVJiKSC1qqC6ISeM3dtwKjgLmhfS4wOkyPAuZ54kWgo5l1baH6RERaXEsF8PXAI2G6i7vvDNO7gC5huhuwLWOb7aFNRKQgpR7AZnYqcBXwaMNl7u6An+T+JppZtZlV19bWNlOVIiItryWOgK8EVrn77jC/u75rITzvCe07gB4Z23UPbUdx99nuXuHuFcXFxSmWLSKSrpYI4C/zUfcDwFJgQpieACzJaB8fzoYYDNRldFWIiBScVO8JZ2anAUOBWzKapwOLzOxmYCtwXWhfBgwHNpOcMXFTmrWJiMSWagC7+36gc4O2t0jOimi4rgO3pVmPiEgu0ZVwIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikaR6R4xcVTL1yVT2u2X6iFT2KyKFSUfAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEkmqAWxmHc2sysw2mtkGM7vIzM40s+Vmtik8dwrrmpnNNLPNZvaKmQ1IszYRkdjSPgKeAfza3c8DyoANwFRghbuXAivCPMCVQGl4TARmpVybiEhUqQWwmZ0BfBF4AMDdD7r7O8AoYG5YbS4wOkyPAuZ54kWgo5l1Tas+EZHY0jwC7gXUAv9lZn8ws/vN7DSgi7vvDOvsArqE6W7Atoztt4e2o5jZRDOrNrPq2traFMsXEUlXmgHcBhgAzHL3C4H9fNTdAIC7O+Ans1N3n+3uFe5eUVxc3GzFioi0tDQDeDuw3d1fCvNVJIG8u75rITzvCct3AD0ytu8e2kREClJqAezuu4BtZtY7NFUCrwJLgQmhbQKwJEwvBcaHsyEGA3UZXRUiIgUn7ZtyfgOYb2anAn8GbiIJ/UVmdjOwFbgurLsMGA5sBg6EdUVEClaqAezuq4GKYyyqPMa6DtyWZj0iIrlEV8KJiESiABYRiUQBLCISSdo/wklKSqY+mdq+t0wfkdq+ReQjOgIWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEkmqAWxmW8xsrZmtNrPq0HammS03s03huVNoNzObaWabzewVMxuQZm0iIrG1xBHwZe7e390rwvxUYIW7lwIrwjzAlUBpeEwEZrVAbSIi0cToghgFzA3Tc4HRGe3zPPEi0NHMukaoT0SkRaQdwA78PzOrMbOJoa2Lu+8M07uALmG6G7AtY9vtoU1EpCC1SXn//+juO8zsbGC5mW3MXOjubmZ+MjsMQT4R4Jxzzmm+SkVEWliqR8DuviM87wEWA4OA3fVdC+F5T1h9B9AjY/Puoa3hPme7e4W7VxQXF6dZvohIqlILYDM7zcw+VT8NfAlYBywFJoTVJgBLwvRSYHw4G2IwUJfRVSEiUnDS7ILoAiw2s/rX+YW7/9rMfg8sMrObga3AdWH9ZcBwYDNwALgpxdpERKJLLYDd/c9A2THa3wIqj9HuwG1p1SMikmt0JZyISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEklWAWxmfdMuRESktcl2QPb/NLNPAA8C8929Lr2SJFeUTH0ytX1vmT4itX2L5IusjoDd/RLgBpKbZtaY2S/MbGiqlYmIFLis+4DdfRMwDZgCXArMNLONZnZNWsWJiBSybPuA+5nZvcAG4HLgn93982H63hTrExEpWNn2Af8HcD9wh7u/V9/o7m+Y2bRUKhMRKXDZBvAI4D13/wDAzE4B2rn7AXd/KLXqREQKWLZ9wE8Bn8yYbx/aRESkibIN4Hbuvq9+Jky3T6ckEZHWIdsA3m9mA+pnzKwceO8E64uISCOy7QO+HXjUzN4ADPg0MC6tokREWoOsAtjdf29m5wG9Q9Mf3f1QemWJiBS+bI+AAQYCJWGbAWaGu89LpSoRkVYgqwA2s4eAzwKrgQ9CswMKYBGRJsr2CLgCON/dPc1iRERak2zPglhH8sObiIg0k2yPgM8CXjWzl4H36xvd/arGNjSzIqAa2OHuI82sF7AA6AzUADe6+8Ew3OU8oBx4Cxjn7ltO5s2IiOSTbAP4ro/xGpNIBvE5Pcz/CLjX3ReY2U+Bm4FZ4fltd/+cmV0f1tOpbiJSsLIdD/hZYAvQNkz/HljV2HZm1p1kHIn7w7yRjKBWFVaZC4wO06PCPGF5ZVhfRKQgZTsc5b+QhOLPQlM34PEsNv0J8F3gwzDfGXjH3Q+H+e1hX/X73AYQlteF9RvWMtHMqs2sura2NpvyRURyUrY/wt0GXAy8C0cGZz/7RBuY2Uhgj7vXfKwKG3D32e5e4e4VxcXFzblrEZEWlW0f8PvhhzIAzKwNyXnAJ3IxcJWZDQfakfQBzwA6mlmbcJTbHdgR1t9Bcsuj7WH/Z5D8GCciUpCyPQJ+1szuAD4Z7gX3KPB/T7SBu3/P3bu7ewlwPfC0u98ArATGhNUmAEvC9NIwT1j+tM47FpFClm0ATwVqgbXALcAykvvDNcUU4Ntmtpmkj/eB0P4A0Dm0fzu8pohIwcp2MJ4PgZ+Hx0lz92eAZ8L0n4FBx1jnb8DYpuxfRCQfZTsWxOsco8/X3f+h2SsSEWklTmYsiHrtSI5Uz2z+ckREWo9sL8R4K+Oxw91/QnKBhYiINFG2XRADMmZPITkiPpmxhEVEpIFsQ/R/Z0wfJrks+bpmr0ZEpBXJ9iyIy9IuRESktcm2C+LbJ1ru7j9unnJERFqPkzkLYiDJ1WoA/wy8DGxKoygRkdYg2wDuDgxw978CmNldwJPu/pW0ChMRKXTZXorcBTiYMX8wtImISBNlewQ8D3jZzBaH+dF8NHi6iIg0QbZnQXzfzH4FXBKabnL3P6RXlohI4cu2CwKgPfCuu88gGbO3V0o1iYi0CtnekujfSIaR/F5oags8nFZRIiKtQbZHwFcDVwH7Adz9DeBTaRUlItIaZBvAB8PdKRzAzE5LryQRkdYh2wBeZGY/I7mf278AT9HEwdlFRCTR6FkQltyJcyFwHsldkXsD/+ruy1OuTUSkoDUawO7uZrbM3fsCCl0RkWaSbRfEKjMbmGolIiKtTLZXwn0B+IqZbSE5E8JIDo77pVWYiEihO2EAm9k57v4X4IoWqkdEpNVo7Aj4cZJR0Laa2WPufm0L1CQi0io01gdsGdO6Bb2ISDNqLID9ONMiIvIxNdYFUWZm75IcCX8yTMNHP8Kdnmp1IiIF7IQB7O5FLVWIiEhrczLDUYqISDNSAIuIRKIAFhGJJLUANrN2Zvayma0xs/Vmdndo72VmL5nZZjNbaGanhvZPhPnNYXlJWrWJiOSCNI+A3wcud/cyoD8wzMwGAz8C7nX3zwFvAzeH9W8G3g7t94b1REQKVmoB7Il9YbZteDhwOVAV2ueS3GEZYBQf3Wm5CqgMQ2GKiBSkVPuAzazIzFYDe0iGsnwNeMfdD4dVtgPdwnQ3YBtAWF4HdD7GPieaWbWZVdfW1qZZvohIqlINYHf/wN37A92BQSSDun/cfc529wp3ryguLv64uxMRiaZFzoJw93eAlcBFJLc1qr8ApDuwI0zvAHoAhOVnAG+1RH0iIjGkeRZEsZl1DNOfBIYCG0iCeExYbQKwJEwvDfOE5U+HG4GKiBSkbAdkb4quwFwzKyIJ+kXu/oSZvQosMLN7gD8AD4T1HwAeMrPNwF7g+hRrExGJLrUAdvdXgAuP0f5nkv7ghu1/A8amVY+ISK7RlXAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkaV6IIdKokqlPprbvLdNHpLZvkeagI2ARkUgUwCIikSiARRpRUlKCmRXso6SkJPZH3GqpD1ikEVu3bqWQB+bTjWfi0RGwiEgkCmARkUgUwCIikSiARZpo+/btjBo1itLSUj772c8yadIkDh48yIMPPsjXv/712OX9nQ4dOsQuQRpQAIs0gbtzzTXXMHr0aDZt2sSf/vQn9u3bx5133pnK6x0+fLjxlSTvKIBFmuDpp5+mXbt23HTTTQAUFRVx7733MmfOHA4cOMC2bdsYMmQIpaWl3H333QDs37+fESNGUFZWxgUXXMDChQsBqKmp4dJLL6W8vJwrrriCnTt3AjBkyBBuv/12Kioq+P73v0/Pnj358MMPj+yrR48eHDp0iNdee41hw4ZRXl7OJZdcwsaNGwF4/fXXueiii+jbty/Tpk1r6Y9IsqDT0ESaYP369ZSXlx/Vdvrpp3POOedw+PBhXn75ZdatW0f79u0ZOHAgI0aMYOvWrXzmM5/hySeTy6/r6uo4dOgQ3/jGN1iyZAnFxcUsXLiQO++8kzlz5gBw8OBBqqurAVi1ahXPPvssl112GU888QRXXHEFbdu2ZeLEifz0pz+ltLSUl156iVtvvZWnn36aSZMm8bWvfY3x48dz3333tewHJFlRAIukYOjQoXTu3BmAa665ht/97ncMHz6cyZMnM2XKFEaOHMkll1zCunXrWLduHUOHDgXggw8+oGvXrkf2M27cuKOmFy5cyGWXXcaCBQu49dZb2bdvH88//zxjx350O8X3338fgOeee47HHnsMgBtvvJEpU6ak/r7l5CiARZrg/PPPp6qq6qi2d999l7/85S+0adPm7y5uMDPOPfdcVq1axbJly5g2bRqVlZVcffXV9OnThxdeeOGYr3Paaacdmb7qqqu444472Lt3LzU1NVx++eXs37+fjh07snr16mNur4sscpv6gEWaoLKykgMHDjBv3jwgOXKdPHkyX/3qV2nfvj3Lly9n7969vPfeezz++ONcfPHFvPHGG7Rv356vfOUrfOc732HVqlX07t2b2traIwF86NAh1q9ff8zX7NChAwMHDmTSpEmMHDmSoqIiTj/9dHr16sWjjz4KJD8OrlmzBoCLL76YBQsWADB//vy0PxJpAgWwSBOYGYsXL+bRRx+ltLSUc889l3bt2vGDH/wAgEGDBnHttdfSr18/rr32WioqKli7di2DBg2if//+3H333UybNo1TTz2VqqoqpkyZQllZGf379+f5558/7uuOGzeOhx9++Kiuifnz5/PAAw9QVlZGnz59WLJkCQAzZszgvvvuo2/fvuzYsSPdD0SaxPL5GveKigqv/4HiZKQ1Bm3m+LNpv0ZLjKNbKK/RmPp/ph/vb8HMCn4siEJ+f7nAzGrcvaJhu46ARUQiUQCLiESiABYRiUQBLNLChg8fzjvvvBO7DMkBOg9YpIW4O+7OsmXLYpciOUJHwCInaerUqUdd2nvXXXdxzz33UFlZyYABA+jbt++RU8G2bNlC7969GT9+PBdccAHbtm2jpKSEN998E4DRo0dTXl5Onz59mD179pF9dujQgTvvvJOysjIGDx7M7t27Adi9ezdXX301ZWVllJWVHTll7eGHHz5yitstt9zCBx980FIfh3wMCmCRkzRu3DgWLVp0ZH7RokVMmDCBxYsXs2rVKlauXMnkyZOPnNq1adMmbr31VtavX0/Pnj2P2tecOXOoqamhurqamTNn8tZbbwHJYDuDBw9mzZo1fPGLX+TnP/85AN/85je59NJLWbNmDatWraJPnz5s2LCBhQsX8txzz7F69WqKiop04UWeSK0Lwsx6APOALoADs919hpmdCSwESoAtwHXu/rYlJ2POAIYDB4CvuvuqtOoTaaoLL7yQPXv28MYbb1BbW0unTp349Kc/zbe+9S1++9vfcsopp7Bjx44jR609e/Zk8ODBx9zXzJkzWbx4MQDbtm1j06ZNdO7cmVNPPZWRI0cCUF5ezvLly4FkFLb6q++Kioo444wzeOihh6ipqWHgwIEAvPfee5x99tmpfgbSPNLsAz4MTHb3VWb2KaDGzJYDXwVWuPt0M5sKTAWmAFcCpeHxBWBWeBbJOWPHjqWqqopdu3Yxbtw45s+fT21tLTU1NbRt25aSkhL+9re/AUeP55DpmWee4amnnuKFF16gffv2DBky5Mg2bdu2PXKBSFFR0QnHA3Z3JkyYwA9/+MNmfpeSttS6INx9Z/0RrLv/FdgAdANGAXPDanOB0WF6FDDPEy8CHc2sKyI5aNy4cSxYsICqqirGjh1LXV0dZ599Nm3btmXlypVs3bq10X3U1dXRqVMn2rdvz8aNG3nxxRcb3aayspJZs2YByfgTdXV1VFZWUlVVxZ49ewDYu3dvVq8v8bVIH7CZlQAXAi8BXdx9Z1i0i6SLApJw3pax2fbQ1nBfE82s2syqa2tr0yta5AT69OnDX//6V7p160bXrl254YYbqK6upm/fvsybN4/zzjuv0X0MGzaMw4cP8/nPf56pU6cet5si04wZM1i5ciV9+/alvLycV199lfPPP5977rmHL33pS/Tr14+hQ4ceGdRdclvqp6GZWQfgMeB2d383c3g8d3czO6mL0N19NjAbkrEgmrNWkZOxdu3aI9NnnXXWcYeUXLdu3VHzW7ZsOTL9q1/96pjb7Nu378j0mDFjGDNmDABdunQ5coZFpnHjxh01QI/kh1SPgM2sLUn4znf3X4bm3fVdC+F5T2jfAfTI2Lx7aBMRKUipBXA4q+EBYIO7/zhj0VJgQpieACzJaB9vicFAXUZXhYhIwUmzC+Ji4EZgrZmtDm13ANOBRWZ2M7AVuC4sW0ZyCtpmktPQbkqxNhGR6FILYHf/HXC8+6FUHmN9B25Lqx4RkVyjsSBEGtGzZ8+Cvrdaw6vzpOUogEUakXnWgkhz0lgQIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEolGQ5OCVzL1yWZdL9OW6SNOehuRejoCFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikaQWwGY2x8z2mNm6jLYzzWy5mW0Kz51Cu5nZTDPbbGavmNmAtOoSEckVaR4BPwgMa9A2FVjh7qXAijAPcCVQGh4TgVkp1iUikhNSC2B3/y2wt0HzKGBumJ4LjM5on+eJF4GOZtY1rdpERHJBS98Ro4u77wzTu4AuYbobsC1jve2hbScNmNlEkqNkzjnnnPQqFTkJTbmbRjZ0x43CFu1HOHd3wJuw3Wx3r3D3iuLi4hQqExFpGS0dwLvruxbC857QvgPokbFe99AmIlKwWjqAlwITwvQEYElG+/hwNsRgoC6jq0JEpCCl1gdsZo8AQ4CzzGw78G/AdGCRmd0MbAWuC6svA4YDm4EDwE1p1SUikitSC2B3//JxFlUeY10HbkurFhGRXKQr4UREIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRtPRdkUWkiXTn5cKjI2ARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCS6EENEjtDFHi1LR8AiIpEogEVEIlEAi4hEogAWEYlEASwiEklOnQVhZsOAGUARcL+7T49ckog0M51p8ZGcOQI2syLgPuBK4Hzgy2Z2ftyqRETSk0tHwIOAze7+ZwAzWwCMAl6NWpWI5J18Oco2d2/WHTaVmY0Bhrn7/wjzNwJfcPevN1hvIjAxzPYG/ngSL3MW8GYzlBtTvr+HfK8f9B5yQb7V39Pdixs25tIRcFbcfTYwuynbmlm1u1c0c0ktKt/fQ77XD3oPuSDf66+XM33AwA6gR8Z899AmIlKQcimAfw+UmlkvMzsVuB5YGrkmEZHU5EwXhLsfNrOvA78hOQ1tjruvb+aXaVLXRY7J9/eQ7/WD3kMuyPf6gRz6EU5EpLXJpS4IEZFWRQEsIhJJqwlgMxtmZn80s81mNjV2PY0xsx5mttLMXjWz9WY2KbSfaWbLzWxTeO4Uu9bGmFmRmf3BzJ4I873M7KXwXSwMP7rmJDPraGZVZrbRzDaY2UX59h2Y2bfCf0PrzOwRM2uX69+Bmc0xsz1mti6j7ZifuyVmhvfyipkNiFf5yWkVAZynlzkfBia7+/nAYOC2UPNUYIW7lwIrwnyumwRsyJj/EXCvu38OeBu4OUpV2ZkB/NrdzwPKSN5H3nwHZtYN+CZQ4e4XkPzAfT25/x08CAxr0Ha8z/1KoDQ8JgKzWqjGj8/dC/4BXAT8JmP+e8D3Ytd1ku9hCTCU5Mq/rqGtK/DH2LU1Und3kj+Wy4EnACO5gqnNsb6bXHoAZwCvE36szmjPm+8A6AZsA84kOevpCeCKfPgOgBJgXWOfO/Az4MvHWi/XH63iCJiP/iOstz205QUzKwEuBF4Curj7zrBoF9AlVl1Z+gnwXeDDMN8ZeMfdD4f5XP4uegG1wH+FLpT7zew08ug7cPcdwL8DfwF2AnVADfnzHWQ63ueet3/frSWA85aZdQAeA25393czl3nyv/ucPY/QzEYCe9y9JnYtTdQGGADMcvcLgf006G7Ig++gE8mgVr2AzwCn8ff/tM87uf65Z6u1BHBeXuZsZm1Jwne+u/8yNO82s65heVdgT6z6snAxcJWZbQEWkHRDzAA6mln9RUC5/F1sB7a7+0thvookkPPpO/gn4HV3r3X3Q8AvSb6XfPkOMh3vc8/Lv29oPQGcd5c5m5kBDwAb3P3HGYuWAhPC9ASSvuGc5O7fc/fu7l5C8pk/7e43ACuBMWG1nH0P7r4L2GZmvUNTJcnwqHnzHZB0PQw2s/bhv6n695AX30EDx/vclwLjw9kQg4G6jK6K3Ba7E7qlHsBw4E/Aa8CdsevJot5/JPkn1ivA6vAYTtKHugLYBDwFnBm71izfzxDgiTD9D8DLwGbgUeATses7Qd39gerwPTwOdMq37wC4G9gIrAMeAj6R698B8AhJn/Uhkn+J3Hy8z53kh937wt/2WpIzPqK/h2weuhRZRCSS1tIFISKScxTAIiKRKIBFRCJRAIuIRKIAFhGJRAEseSuMFndFg7bbzSyrwVjM7H+Z2T+lU51I43QamuQtM5sIXOTuN2W0vQh8191/28i2Re7+Qdo1ipyIjoAln1UBI+rHsg2DFn2GZLjR6jAG7t31K5vZFjP7kZmtAsaa2YNmNiYs+1cz+30YM3d2uGoMM3smbPOymf3JzC4J7UVm9u9h/VfM7BuhvdzMnjWzGjP7Tf2lsyLHogCWvOXue0mu5royNF0PLCK50rEC6Adcamb9MjZ7y90HuPuCBrv7P+4+0JMxcz8JjMxY1sbdBwG3A/8W2iaSDJfY3937AfPD2B3/AYxx93JgDvD95nm3UogUwJLvHiEJXsLzI8B14Sj3D0AfkkH46y08zn4uC3eIWEsyaFCfjGX1AyHVkIQuJIPc/MzDkI7hfwa9gQuA5Wa2GphGMjCMyDHlzG3pRZpoCXBvuA1Ne2Av8D+Bge7+tpk9CLTLWH9/wx2YWTvgP0nGENhmZnc12Ob98PwBJ/6bMWC9u1/UxPcirYyOgCWvufs+kpG95pAc/Z5OErJ1ZtaFj7onTqQ+bN8M4y+POdHKwXLglvohHc3sTJI7MRSb2UWhra2Z9TnBPqSVUwBLIXiE5H5tj7j7GpKuh43AL4DnGtvY3d8Bfk4yWthvSIYvbcz9JEM9vmJma4D/7u4HScL7R6FtNfDfTvbNSOuh09BERCLREbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhLJ/wfCfPvQ3t9jHwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "random.seed(1)\n",
    "perm_variance = [perm_test(four_sessions) for _ in range(3000)]\n",
    "print('Pr(Prob)', np.mean([var > observed_variance for var in perm_variance]))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.hist(perm_variance, bins=11, rwidth=0.9)\n",
    "ax.axvline(x = observed_variance, color='black', lw=2)\n",
    "ax.text(60, 200, 'Observed\\nvariance', bbox={'facecolor':'white'})\n",
    "ax.set_xlabel('Variance')\n",
    "ax.set_ylabel('Frequency')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## F-Statistic\n",
    "We can compute an ANOVA table using statsmodel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            df  sum_sq     mean_sq         F    PR(>F)\n",
      "Page       3.0   831.4  277.133333  2.739825  0.077586\n",
      "Residual  16.0  1618.4  101.150000       NaN       NaN\n"
     ]
    }
   ],
   "source": [
    "model = smf.ols('Time ~ Page', data=four_sessions).fit()\n",
    "                \n",
    "aov_table = sm.stats.anova_lm(model)\n",
    "print(aov_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F-Statistic: 1.3699\n",
      "p-value: 0.0388\n"
     ]
    }
   ],
   "source": [
    "res = stats.f_oneway(four_sessions[four_sessions.Page == 'Page 1'].Time, \n",
    "                     four_sessions[four_sessions.Page == 'Page 2'].Time,\n",
    "                     four_sessions[four_sessions.Page == 'Page 3'].Time,\n",
    "                     four_sessions[four_sessions.Page == 'Page 4'].Time)\n",
    "print(f'F-Statistic: {res.statistic / 2:.4f}')\n",
    "print(f'p-value: {res.pvalue / 2:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two-way anova only available with statsmodels\n",
    "```\n",
    "formula = 'len ~ C(supp) + C(dose) + C(supp):C(dose)'\n",
    "model = ols(formula, data).fit()\n",
    "aov_table = anova_lm(model, typ=2)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chi-Square Test\n",
    "## Chi-Square Test: A Resampling Approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Headline  Headline A  Headline B  Headline C\n",
      "Click                                       \n",
      "Click             14           8          12\n",
      "No-click         986         992         988\n"
     ]
    }
   ],
   "source": [
    "# Table 3-4\n",
    "click_rate = pd.read_csv(CLICK_RATE_CSV)\n",
    "clicks = click_rate.pivot(index='Click', columns='Headline', values='Rate')\n",
    "print(clicks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Headline A</th>\n",
       "      <th>Headline B</th>\n",
       "      <th>Headline C</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Click</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Click</th>\n",
       "      <td>11.333333</td>\n",
       "      <td>11.333333</td>\n",
       "      <td>11.333333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>No-click</th>\n",
       "      <td>988.666667</td>\n",
       "      <td>988.666667</td>\n",
       "      <td>988.666667</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          Headline A  Headline B  Headline C\n",
       "Click                                       \n",
       "Click      11.333333   11.333333   11.333333\n",
       "No-click  988.666667  988.666667  988.666667"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Table 3-5\n",
    "row_average = clicks.mean(axis=1)\n",
    "pd.DataFrame({\n",
    "    'Headline A': row_average,\n",
    "    'Headline B': row_average,\n",
    "    'Headline C': row_average,\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed chi2: 1.6659\n",
      "Resampled p-value: 0.4820\n"
     ]
    }
   ],
   "source": [
    "# Resampling approach\n",
    "box = [1] * 34\n",
    "box.extend([0] * 2966)\n",
    "random.shuffle(box)\n",
    "\n",
    "def chi2(observed, expected):\n",
    "    pearson_residuals = []\n",
    "    for row, expect in zip(observed, expected):\n",
    "        pearson_residuals.append([(observe - expect) ** 2 / expect\n",
    "                                  for observe in row])\n",
    "    # return sum of squares\n",
    "    return np.sum(pearson_residuals)\n",
    "\n",
    "expected_clicks = 34 / 3\n",
    "expected_noclicks = 1000 - expected_clicks\n",
    "expected = [34 / 3, 1000 - 34 / 3]\n",
    "chi2observed = chi2(clicks.values, expected)\n",
    "\n",
    "def perm_fun(box):\n",
    "    sample_clicks = [sum(random.sample(box, 1000)),\n",
    "                     sum(random.sample(box, 1000)),\n",
    "                     sum(random.sample(box, 1000))]\n",
    "    sample_noclicks = [1000 - n for n in sample_clicks]\n",
    "    return chi2([sample_clicks, sample_noclicks], expected)\n",
    "\n",
    "perm_chi2 = [perm_fun(box) for _ in range(2000)]\n",
    "\n",
    "resampled_p_value = sum(perm_chi2 > chi2observed) / len(perm_chi2)\n",
    "print(f'Observed chi2: {chi2observed:.4f}')\n",
    "print(f'Resampled p-value: {resampled_p_value:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed chi2: 1.6659\n",
      "p-value: 0.4348\n"
     ]
    }
   ],
   "source": [
    "chisq, pvalue, df, expected = stats.chi2_contingency(clicks)\n",
    "print(f'Observed chi2: {chi2observed:.4f}')\n",
    "print(f'p-value: {pvalue:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure chi-sq distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR8AAACsCAYAAABCfqgiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2oUlEQVR4nO2deXxMV//HPyf7HmQVIQQhiRIkVYJaa6laSh9Ka+9GVUu1lqYtpbTFr2p5KlSpLYpaY6m1VFFBbEFJBEmzIrLIPp/fHzOZJ5Ftkkwykzjv1+u+cu+5Z/nOzdzPnPV7BElIJBJJVWOgawMkEsmziRQfiUSiE6T4SCQSnSDFRyKR6AQpPhKJRCcY6dqAp7G3t2fDhg11bYZEItEC58+fTyTpUNQ9vROfhg0bIiQkRNdmSCQSLSCEuFvcPdnskkgkOkGKj0Qi0QlSfCQSiU7QqM9HCNEbwBIAhgBWk1zw1P13AUwEkAsgFcDbJMNU92YAGKe69wHJg9ozXyKpGWRnZyMqKgoZGRm6NqVcmJmZwdXVFcbGxhqnKVV8hBCGAJYD6AkgCsA5IcTuPHFRsYnkj6r4/QEsBtBbCOEFYBgAbwAuAA4LITxI5mpsYTFkZWXh5s2bcHFxgZ2dXUWzk0h0SlRUFKytrdGwYUMIIXRtTpkgiQcPHiAqKgqNGjXSOJ0mza7nAdwmGUEyC0AQgAFPFZ6c79ISQN5q1QEAgkhmkrwD4LYqvwoTERGBli1bYv/+/drITiLRKRkZGbCzs6t2wgMAQgjY2dmVudamSbOrHoD7+a6jALQrwoCJAKYAMAHQLV/aM0+lrVdE2rcBvA0ADRo00MRuNGrUCEFBQfD399covkSi71RH4cmjPLZrrcOZ5HKSjQF8CuCzMqYNJOlL0tfBocj5SIUwNTWFn58fHj9+XA5rJRKJrtGk5hMNoH6+a1dVWHEEAfhvOdOWiTfffBOmpqY4evSotrKUSGokX375JaysrNCvXz8MGzYMQghs27YNjRs31plNmojPOQBNhRCNoBSOYQCG548ghGhK8pbq8mUAeee7AWwSQiyGssO5KYC/tWE4ACxYsABmZmbayk4iqfHs3LkTQ4YMwWeflalxUimUKj4kc4QQ7wM4COVQ+xqS14QQcwCEkNwN4H0hRA8A2QAeARilSntNCPErgDAAOQAmamOkKw8vLy9cvnwZOTk5MDLSu5UiEolOmTdvHtatWwdHR0fUr18fnp6eWLFiBQwNDXHkyBEcO3ZMp/Zp9MaS3Adg31Nhn+c7n1xC2nkA5pXXwJLYuXMnxo8fj4iIiDIN8UkkNZ3z588jKCgIoaGhyMnJQZs2bdC2bVu8++67sLKywscff6xrE/VvYWlZ6N27Nw4fPgwnJyddmyKR6BUnT57EoEGDYGFhAQDo37+/ji0qTLVeXuHk5ARzc3PExMTo2hSJRFJGqrX4kESnTp2wdu1aXZsikegVnTt3xs6dO5Geno6UlBTs2bNH1yYVolo3u4yNjXHo0CF4eHjo2hSJRK9o06YNhg4dilatWsHR0RF+fn66NqkQ1Vp8AMDFxQW3bt2Cq6urrk2RSPSKWbNmYdasWbo2o1iqdbMLAL755hu8+eabujZDIpGUkWovPtOnT8e+fftKjyiRSPSKai8+Dg4OiIyMxIMHD3RtikQiKQPVXnyuXbuGAQMG4Pz587o2RSKRlIFq3+HcunVr/P333/Dy8tK1KRKJpAxU+5qPlZUVUlJScOnSJV2bIpFIykC1Fx8AeP/99/Hdd9/p2gyJpFozduxYODo6okWLFlVSnkbiI4ToLYS4KYS4LYSYXsT9KUKIMCHEZSHEESGEW757uUKIUNWxW5vG57Ft2zasXLmyMrKWSJ4ZRo8ejQMHDlRZeaWKTz4H8n0AeAF4XeUYPj8XAfiSbAlgG4Bv891LJ+mjOipldZupqSl2796N3FyteeuQSJ45OnfujDp16lRZeZp0OKsdyAOAECLPgbx69wqS+R2DnAHwhjaNLI1jx47hrbfeQvfu3aVrDUm158MPP0RoaKhW8/Tx8cH333+v1TwriibNrqIcyBdyAp+PcQDybylhJoQIEUKcEUIMLCqBEOJtVZyQhIQEDUwqyODBg3Hnzh24ubmVHlkikegFWh1qF0K8AcAXwIv5gt1IRgsh3AEcFUJcIRmePx3JQACBAODr60uUERsbGwQHB6Np06Zo167QxhoSSbVC32oolYUmNR+NnMCr3KjOAtCfZGZeOMlo1d8IAMcBtK6AvYUgCQMDA0yYMAEbNmzQZtYSiaQS0UR81A7khRAmUDqQLzBqJYRoDWAllMITny+8thDCVHVuD8Af+fqKKsKpU6dgY2ODEydOQAiBS5cuYeHChdrIWiJ5Jnn99dfRvn173Lx5E66urvjpp58qtTxtOZD/DoAVgK2qzcPuqUa2PAGsFEIooBS6BU9ts1xubGxskJKSgri4OABAfHw89u3bh4kTJ2oje4nkmWPz5s1VWp62HMj3KCbdXwCeq4iBxZHntzlPfPbt24d58+bhrbfegomJSWUUKZFItEi1neFsZ2cHQ0NDtfh89NFHSEpKksIjkVQTqq34GBoawsHBQS0+QgisWrVK6/MjJBJJ5VBtxQdQNr1iY2MBKEe9pkyZovON0CQSiWZUa5caTk5O6ppPrVq1kJiYCDs7Ox1bJZFINKFa13ycnZ3V4gMAf/31F5YsWaJDiyQSiaZUa/HJq/mQyknRwcHBWLx4sY6tkkiqJ/fv30fXrl3h5eUFb2/vSv8hr/bik5mZieTkZADAokWLcOfOHR1bJZFUT4yMjLBo0SKEhYXhzJkzWL58OcLCtDItr0iqvfgAUHc6p6WlYdasWbh48aIuzZJIqiV169ZFmzZtAADW1tbw9PREdHShlVRao1qLj7OzMwAUGG5fuHChdKkqqfZ06dJFvQ24ts7LQmRkJC5evFipC7Wr/WgX8D/xcXBwQFpampxoKJFUgNTUVAwePBjff/89bGxsKq2cGiU+gNKlalhYGObOnasrsySSCnP8+HGtn2tCdnY2Bg8ejBEjRuDVV18tU9qyUq2bXXZ2djAwMCggPn///Td27dqlHgGTSCSaQRLjxo2Dp6cnpkyZUunlVYUD+VFCiFuqY5Q2jTc0NISjo6O6wxkAFi9ejCtXrkC1ul4ikWjIqVOnsH79ehw9ehQ+Pj7w8fGp1K3IS2125XMg3xNKF6rnhBC7n3KNkedA/okQ4j0oHcgPFULUAfAFlN4NCeC8Ku0jbX2A/LOcAaVrjc8++wxjxoyBv7+/toqRSGo8HTt2rNIWgyY1H7UDeZJZAPIcyKsheYzkE9XlGSi9HQJALwCHSD5UCc4hAL21Y7qSp8XH3Nwcu3fvRmRkpDaLkUgkWkaTDueiHMiXNP6W34F8WZ3PlxknJyfcvHlTfW1ra4v4+PgSUkgkEn1Aqx3O+RzIl2n70IrsXuHs7IzY2NgC1cXAwEC88847ZcpHIpFULZXtQF6jtCQDSfqS9HVwcNDUdgCFl1gAQFRUFG7cuCFHvCQSPUaTZpfagTyUwjEMwPD8EfI5kO+d34E8lH6fvxZC1FZdvwRgRoWtzkf+uT62trYAgDlz5mizCIlEUgmUWvMhmQMgz4H8dQC/5jmQF0LkbX+c34G8ek92kg8BfAWlgJ0DMEcVpjWKmmiYlJSEAQMGYNu2bdosSiKRaJFKdSCvurcGwJryGlgaT6/vApQ7W0RHRyM1NbWyipVIaiQNGzaEtbU1DA0NYWRkhJCQkEorq1ovrwAKr2wHAAMDg0p9aBJJTebYsWOwt7ev9HKq9fIK4H+7WMTExBQI37RpE9q2bYucnBwdWSaRSEqi2ouPoaEhGjRoUGhSoYWFBVxcXJCUlKQTuySSitClSxf1AQBr165VX+d3l6FpHE0RQuCll15C27ZtERgYqL0PVATVvtkFAI0aNUJERESBsIEDB2LgwIG6MUgiqab8+eefqFevHuLj49GzZ080b94cnTt3rpSyhL7NhfH19WVZ+2veeust7Nmzp0C/DwD06tULPj4++Oabb7RpokSida5fvw5PT09dm1GAL7/8ElZWVvj44481il/UZxBCnCfpW1T8at/sApQ1n7i4OKSlpRUIb9q0KVxdXYtJJZFI8pOWloaUlBT1+e+//44WLVpUWnk1otnl7u4OQOn60dvbWx2+bNkykARJ6WJDIimFuLg4DBo0CACQk5OD4cOHo3dvra4DL0CNEJ9GjRoBAO7cuVNAfA4fPowRI0bg6NGjBcIlEklh3N3dq9T/eY1pdgEo1OncqFEj9O3bF0ZGNUJjJZIaRY14Kx0cHGBpaVloz67GjRvj559/RkZGho4sk0gkxVEjaj5CiCKH2wFgwoQJ8PHxqXqjJBJJidSImg+gbGIVtVtpjx490KBBAygUChgY1AitlUhqBDVGfNzd3XHs2LFCI1uvvvqq2t9PrVq1dGegRCIpgLZ2r+gshLgghMgRQgx56l6uys2G2tVGZdCoUSOkpqYiMTGxQHh2djYcHR0xf/58AMCFCxcAAGvWrMGIESMQGBiI+/fvF8pPIpFULqWKT77dK/oA8ALwuhDC66lo9wCMBrCpiCzSSfqojv5F3NcK+Yfb82NsbIw5c+agQ4cOGD16NPz8/JCRkYHs7GwcPXoU77zzDn777TcoFAqU1YWrRFKTGDt2LBwdHQtMLHz48CF69uyJpk2bomfPnnj0SGsbz2ht94pIkpcBKLRmWRnJm2hYVKfzpEmTEBAQgI0bN2L69OkQQuCdd97Bv//+ixs3bmD48OFYu3YtmjRpgl9//bWqTZdI9ILRo0fjwIEDBcIWLFiA7t2749atW+jevTsWLFigtfI0EZ+K7kBhpnIOf0YIMbCoCBVxIJ9Hw4YNARSu+SQnJyMoKAhXrlzBmjVrMG/ePJiamuaVi2bNmsHBwQH+/v5o0aIFhg4dip07d5bLBomkOtO5c2fUqVOnQNiuXbswapRyr89Ro0Zp9d2oig5nN5LRQgh3AEeFEFdIhuePQDIQQCCgXFhankKsrKzg4OBQoOaTm5uLESNGIDw8HCtWrMBLL71UbPpmzZrh6NGjWLhwIXr37o3IyEg0aNBAjpBJdMLTbjBGjx6N0aNHlzu8rHu25xEXF4e6desCUHoNze8xtKJobfeK4iAZrfobAeA4gNZlsK9MuLu7Izz8f7q2atUq7N27FxMnTsTrr7+O27dvl5je1NQUs2bNQkJCAlq3bo3JkyfLHTAkEhVCCK2ukdTK7hXFodq14gnJTCGEPQB/KLdSrhS8vb2xd+9e9bW1tTWGDBmCCRMm4MMPP8SqVavw8OFDmJmZlZiPq6srxo4di8WLF6NOnTqYPXt2ZZkskRRJcTUVbYVripOTE2JiYlC3bl3ExMTA0dGxQvnlRyu7Vwgh/IQQUQBeA7BSCHFNldwTQIgQ4hKAYwAWPLXHu1Z57rnnEB8fj7i4OISEhOC1117D1q1bIYTAu+++i4MHD8LY2LjUfIQQWLhwIWbNmiUdkkmeafr3749169YBANatW4cBAwaUkqIM5Lmc0Jejbdu2LC9HjhwhAG7dupXW1tZ8++23C9wPCQnh8ePHy5RnUlISX375ZR45cqTcdkkkpREWFqZrEzhs2DA6OzvTyMiI9erV4+rVq5mYmMhu3bqxSZMm7N69Ox88eFBs+qI+A4AQFvOu15gZzoCy5gMAP/74I9LS0jB58uQC9z/44APk5OTg7NmzZcr3zp07GDx4MEJCQtC4cWOt2SuR6BObN28uMvzIkSOVUl6NGspxcHCAs7MzUlNT8e6778LLq+BcyJUrV2LPnj1lytPW1hZ79uyBgYEBfv/9d22aK5E809Qo8QGA5s2bIysrC8uXLy90r0mTJti/fz/OnTtXpjzd3d1x8+ZNvPfeezh16pQcAZNItECNE5+oqCiEhoYiOzu70D0DAwNMnDixXLOY7e3tERwcjI4dO+Knn37ShqkSyTNNjRKfuLg43LlzBySLnNNjYmKC0NDQcu9m0bt3b/To0QOTJk3CtWvXSk8gkUiKpUaJz40bN9RuM65cuVJkHCsrK8yYMQNhYWUf8Tc0NMSGDRswZswY1K9fv/QEEomkWGqU+Lz44ou4e/cuDA0Ncfny5SLjCCGwZMmScu/l7uTkhBUrVuDSpUuYOXNmRcyVSJ5paoz43L9/H97e3ggJCUGzZs2Krfk4OTnhwYMHGDlyZIXKO3jwIObPn48tW7ZUKB+JRF+4f/8+unbtCi8vL3h7e2PJkiUAKs+tRo0Rn+DgYISFhcHR0REtW7YstuYDALdv30a/fv0q5ETsiy++QPv27TFjxowiO7clkuqGkZERFi1ahLCwMJw5cwbLly9HWFhYpbnVqDHic+DAAbi7u6N58+Zo3bo1IiMji3UOZmpqirCwMNy7d6/c5RkbGyMoKAgnTpxAbm4ucnJyyp2XRKIP1K1bF23atAGgXBfp6emJ6OjoSnOrUWNmOG/atAmRkZEQQsDf3x8AcOrUqSLXZjVr1gzh4eEVXqHboEEDpKSk4IUXXsArr7yCr776qkL5SSR5PO0iI4/jx49j7dq1WLt2baF7xbnXyEtXFiIjI3Hx4kW0a9eu0txq1Iiaz+nTpxEYGAg3NzcAgK+vL0xNTXHy5Mki4wshsHfvXnh4eBTy+VxWrK2t0aZNG8ybN6/CK4glEn0gNTUVgwcPxvfffw8bG5sC97TqVqO4RV/5DwC9AdwEcBvA9CLudwZwAUAOgCFP3RsF4JbqGFVaWeVZWPr222/TysqKmZmZ6rBOnTrx+eefLzZNaGgoX375Zf7zzz9lLu9pUlJS2KpVK/7yyy8VzkvybKIPC0tJMisriy+99BIXLVqkDvPw8OC///5Lkvz333/p4eFRZNqyLiytVAfyQog6AL4A0A5KX9BfqHz8aA2SCA4ORq9evWBiYqIO79SpEy5cuIC0tLQi07Vq1Qp79+7Vin8SKysrhISEYMSIEdi6datcfiGplpDEuHHj4OnpiSlTpqjDK8utRmU7kO8F4BDJhyQfATgEZS1Ka+Tk5CAgIAATJkwoEN6xY8dSV7Bv2bIFDg4OiIyMrLAdRkZG2L59O/7zn/9g2bJlFc5PIqlqTp06hfXr1+Po0aPw8fGBj48P9u3bh+nTp+PQoUNo2rQpDh8+jOnTC+2eVS406XAuyoF8Ow3z18j5vBDibQBvA8pO3LIQGxuLoUOHFtoQsH379hBC4M8//0S3bt2KTPvCCy9gypQpBWpMFWHIkCHo168fPv74Y3Tr1g3e3t4lxk9OToaNjQ127NiB3377DQkJCXjy5AmMjY1x5MgR7N69GxcuXICHhwf8/PzQpEkTrbqxlEjy07Fjx2Jr7ZXhVkMvOpxJBpL0Jenr4OBQprSzZs2Cp6dnoYdWq1YttGzZsthOZwBwc3PD119/jaioKK00lYQQ+Pnnn/H555/Dw8OjyDj37t3DvHnz0LJlS7i4uCArKwuXLl3CiRMn8PDhQxgaGqo7+U6cOIE5c+ZgxIgR8PDwQNOmTUESN2/eRFJSUoXtlUh0SnGdQfxfh3F7AAfzXc8AMKOYuGuRr8MZwOsAVua7Xgng9ZLKK2uHs5ubG4cMGVLkvYkTJ9LS0pLZ2dnFpv/5558JgOfOnStTuaWxefNmjh8/ngqFQh22YcMGCiEIgJ06deL8+fOZkpJSIM7TpKen8/LlywwMDOQ333xDkuzSpQtNTU35n//8h/v27WNubq5WbZdUPfrS4VwRytrhrIn4GAGIANAIgAmASwC8i4n7tPjUAXAHQG3VcQdAnZLKK4v4JCcn08/Pj8uXLy/yflBQEAHwzJkzxeaRlJTEtWvX8smTJxqXqwlff/01AXDu3Ll89dVXuWvXLiYmJjIgIIAREREVyvv06dOcNGkS7ezsCIDjx48nyQKjfZLqRVhYWIk/QvqOQqHQvvgo06MvgH8AhAOYpQqbA6C/6twPyv6cNAAPAFzLl3YslEP0twGMKa2sivhwfprExEQaGBjws88+KzHe48ePuWzZMj58+FBrZWdlZbFZs2YUQtDCwoJr1qzRWt55ZGZmcsuWLTxz5gzv3r1LBwcHBgQEMDExUetlSSqXiIgIJiQkVEsBUigUTEhIKPJHtSTxEdSzYWFfX19quuJ8xowZePDgAQIDA4uN8+KLLyIpKQmXLl0qNk5oaChat26NdevWVXjBaR6TJk3CsmXL0KdPH3z11Vdwc3ODvb29VvIuivDwcEybNg07duyAtbU1pkyZgmnTpsHS0rLSypRoj+zsbERFRSEjI0PXppQLMzMzuLq6FtodRghxnqRvkYmKUyVdHWWp+Xh4ePCVV14pMc7ChQsJgHfu3Ckx3tWrVzUutySOHz/O6Oho3rx5k2vWrGFGRgYbN27MHj16lNj3pC2uXLnCwYMH08XFhWlpabxx44Zsjkl0Bira7KrKQ1PxiYmJIQB+++23Jcb7559/CIA//PBDifESExP53nvv8dixYxqVXxSrVq2ikZERR4wYUSB89erVBMCPP/643HmXlUePHjEzM5Nubm5s0qQJd+3aVS2r9JLqTY0Un9zcXF6+fFk97bskPD092aNHjxLjpKen083NjUuXLtWo/PwoFArOmjWLANi7d28mJSUVijN9+nQeOHCgzHlXBIVCwX379tHT05MA2KNHD60sJ5FINKVGik9Z+PTTT2lkZMRHjx6VGC+veVLWka8TJ04QAN96660Sm1aZmZl85513Shx9qwyysrK4ZMkSOjo6MiwsjNHR0UxOTq5SGyTPJs+8+Jw6dYoAuGnTplLjDho0iAMHDtQo39zcXB4+fJgkefjw4VKbNYmJiXR3d6eTkxMjIyM1KkObpKenkyT79OlDFxcXbtmyRTbFJJVKSeKjFzOcK5t27drB1dUVv/zyS6lxu3Tpgi5duiiVuQRIYvLkyejRowf+/PNPdO/evdSlD3Z2dti7dy8yMjI0skXbmJmZAQC+/PJLODs7Y+jQoejVqxdu3bpV5bZIJDqv6Tx9VEbNhyQDAgIohOC9e/dKjZuVlVVq/8ynn35KAJw6dWqZaw8RERFUKBQ8ceIEMzIyypRWW+Tk5HDp0qW0sbHhli1bmJ6erjNbJDUXPOvNLlL5wgPg7NmzS42bNzx/5cqVIu8HBgYSAN99991yN1tu375NIyMjDhs2TKfLI/Imts2cOZMeHh48cuSIzmwpidzcXD548IA3b95kSEgIjx8/zgMHDvDx48eMjY3lvn37+Pvvv/PPP//kxYsXeefOHaakpOja7GceKT4qunfvzoYNG5b6sicnJ3PPnj1FCktOTg5jYmI4c+ZM5uTkVMieBQsWqJdg6JqDBw+ycePGBMA33niDsbGxVW5DVlYWr169ys2bN/OLL77giBEj1LPTW7VqRQCFjpCQEO7atavIey+//DJJZT9e+/btOWzYME6fPp2rV6/mH3/8ofUlNZLClCQ+NcaHsyaMHz8er7/+Oo4ePYoePXoUG8/a2hrdunXDp59+is6dO6Nfv34AgD///BMTJ07E7t27MW/evArb88knn8DY2BjDhw9HSkoKrK2tK5xneXnppZdw5coVzJ8/HwsWLICbmxu++uorZcegQeV0DUZHR+PkyZPo1KkT7ty5gx49eiAzMxOA0kNA/fr11a5SJk+ejMePH8POzg62trawsrKCmZkZmjVrBnd3d5w+fRrZ2dlIT09HamoqkpKS4OzsDABo1KgRkpKS8Pfff2Pbtm1qZ//Xr19HREQEFi1aBB8fH7Ru3Rq+vr7w8PCotM8syUdxqqSrozJrPunp6axduzYHDx5catysrCx6enryq6++IkmGh4fT3t6eHh4efPDggVbtio+PZ9OmTTlnzhyt5lterl+/ztTUVG7bto1+fn4MCQnRSr5paWncvXs3J0yYQA8PD3UNZdWqVYyPj+fUqVO5YcMGhoaGqkfmtE12djbDw8O5f/9+Zmdnc/fu3fTz86OZmZnaHhsbG3733XckyT/++ENOS6gAkM2u/zFz5kwKITRyYZDXZxAaGkovLy/Wrl27Uibp5ebmcuTIkQSgFjt9YMeOHXR0dKQQghMmTCjXwtsHDx5w8+bNVCgUnDp1KgHQ0tKSffv25aJFi3ju3LkqWXZSGtnZ2bxy5Qp//vlnvvvuu9yxYwcTEhIIgAYGBvT19eXHH3/M4OBgKUZloMLig9IdyJsC2KK6fxZAQ1V4QwDpAEJVx4+llVXZ4hMfH08LCwu+8cYbGsVfvXo1DQ0NaWJiwqNHj1aaXTk5OXzrrbe4Y8cO5ubm6s38m0ePHvGDDz6ggYEBhw4dqlGatLQ0btq0if369aORkREBMDQ0lDdu3OChQ4eqzahaRkYGf//9dwYEBLBz5840MTEhALUD9eDgYP711196IZ76SoXEB4AhlK403PE/fz5eT8WZkCcsAIYB2ML/ic/V0spgFYoPSU6dOpUGBga8detWqXEPHDjA2bNn8/bt25VuF6lcEvHOO+9wwoQJFe7Q1iYXL17krVu3eOPGDXbo0IGnTp0qMt7cuXNpbW1NAHR1deW0adN47tw5vRHTivDkyRMePnyY27dvJ0k2a9aMAFirVi2+9tpr/OmnnzRa7vMsUVHxKdWTIYCDANqrzo0AJAIQ+io+MTExNDMz49ixY0uMl+eMbP369dyyZQsXLFhQ6bYpFAp+8sknBMD//Oc/zMrKqvQyy8KRI0fo4uJCABw+fDjDwsK4fPlydu3alRkZGVy+fDlHjRrFo0eP1ngPi4mJidyyZQvHjh2rfiZCCMbGxvLy5cs8c+ZMjX8GpVFR8RkCYHW+6zcBLHsqzlUArvmuwwHYq8QnDcBFAH8A6FRMGW8DCAEQ0qBBgyp5KJMnT6ahoSFDQ0OLvB8aGkpzc3P6+/szMzOTo0ePpr+/f5WJwXfffccRI0ZQoVDoXbU+JSWF7733Hg0NDWlgYEAAbNOmDcPDw3Vtms5QKBQMDQ3lypUrSVLdh+fo6MjRo0dz+/btz+S8I12KjykAO1VYWyh3srApqbyqqPmQyo5Qe3t7dujQodCvU2JiIhs2bMh69eoxJiaGpLLKnZWVxVu3blXZHBiFQsEjR46wefPmvHbtWpWUqQmfffYZAdDExIS9evXisWPH2LJlSy5dulT6DlKRmJjIDRs28PXXX2etWrXUz+v48eNMSkrSaKZ9TUBnza4i8joOwLek8qpKfEhyzZo1BFDIxWlCQgL79OlTaPV5WloanZ2d2b9//yqz8cyZM3RycqKVlRUPHTpUZeXmJzs7mxs3buQLL7zA2NhYHjx4kF9++aVahO/du8cuXboQAN3d3fnLL7/oVX+VrsnKyuKxY8c4ZcoUJiUlceXKlQRAHx8fBgQE8Ny5czW2eVZR8SnVgTyAiU91OP+qOncAYKg6dwcQDS06kK8oubm59Pf3p729PePj40kqmzvFLasglR3QkZGRTE1NrbIvTFRUFAcMGMCoqCjev3+/ymoXaWlpXLp0Kd3c3AiAnp6exe7ykec7yMfHh8bGxrx79y7j4+Nr7EtVEe7cucNvv/2WHTt2VDdbXVxcGBQURIVCUaNmXldIfJTpS3UgbwZgK5RD7X8DcFeFDwZwDcph9gsAXimtrKoUH1LpdtTU1JS9evXi+vXrCYAfffRRiWmSk5Pp4+PDadOmVZGVSrKzs/ncc8/R19e30p2C7dy5k/b29gTADh06cNeuXRoJSW5uLi9cuEBSuZzFy8uLGzdu1Lt+K30hISGB69at45AhQ3j8+HFevXqV5ubm7NevH1euXMmoqChdm1ghKiw+VXlUtfiQ5I8//kgANDIyYufOnUvtVM4bkTpw4ECVz8n57bffWLt2bTo4ODAtLU2recfGxvLTTz/lpUuXePnyZfbt25cnT54sV14KhYJBQUH08vIiADZu3Jh79uzRqr01kdu3b3PSpEls2LChesZ1q1at1B3Z1W3KghSfUkhPT6eFhQUBcOvWrRqny8nJ4ciRIzl16tRKtK4w9+/f586dO5mbm8upU6fy5s2bFcovMjKSEydOpJmZGYUQpfq7Lgu5ubn87bff2LZtW+7evZtxcXH86quv1M1cSdEoFApeu3aN33zzDV988UXOnTuX2dnZbNy4MYcOHcq1a9eqB0P0GSk+JZCVlUWFQsFdu3axcePGtLa21ngtk0Kh4KRJk9RLIqq6k/Xq1au0tbWliYlJuVbGR0dHc9SoUTQyMqKxsTHHjRtXac05hUJBhULBtWvXEgBNTU05ZswYra0bexZ4+PAhx4wZQ2dn5wK1ooCAAJL6WSuS4lMMubm5fO211/jmm29SoVAwKiqKbm5utLe3540bNzTKI++lWrlyJTt27Fik8/jKJCYmhiNGjGBAQAAVCgX37t1bat/M2bNnefbsWcbFxdHOzo6TJ0+u0qHfsLAwvvfee7SwsGCtWrX45MkTXrt27ZmcB1MeFAoFL168yK+//ppdunRh3759SSpdh/Ts2ZMLFixgSEiIXow4SvEphmnTphGAegUzqdxqx8HBgY6OjmX6Vd6yZQsHDRrE3NxcnfjCyRttAsDWrVtz3759Rd7PGxLv3bs3SVba6nFNSEpK4okTJ6hQKPjcc8/R0tKSY8aM4bFjx+QoWRnIq/HMmDGD3t7e6lpRrVq12L9/f8bFxTEtLU0nYiTFpwh++OEHAuCECRMKVVevX79ONzc3Wlpacv/+/RrnqVAoePnyZVpYWPCXX37Rtsmlkpuby/Xr19Pd3V3t3nXv3r1cvXq1+kvp6urKRYsW6dXKbIVCwb/++otjx46llZUVAXDKlCkkq/8e5rogJiaGGzdu5Pjx49myZUtmZ2dz9uzZtLW1Zd++fTl//nyePHmySn54pPg8hUKhYP/+/TlgwIBih4D//fdf+vj4UAjBgIAAjYeKU1NTOWXKFMbGxjIkJITBwcFV/vJkZWXx9u3baoE1NDSkq6sr161bp/czkNPS0rhhwwaeP3+e169fJwA2adKEn376KU+fPi1rROXk8OHDfPvtt9m8eXN1zcjExISbNm1iamoqf/31V969e1fr31UpPvnYvn07T58+zaysrFKVPzU1lWPGjCEA+vv7a+QDKD+DBg1i/fr1mZGRUSV9QXlO6YcOHUojIyM2atSImzZtYsuWLdU1vNmzZ/Ps2bPVojaRNxv4pZdeUrvmOH78OMPDw7lt27Yq71+rKcTHx3Pnzp2cNm0ar169yqNHj6oFycnJif369ePs2bN5/PjxCpclxUfF+vXraWhoyF69epUp3caNG1mrVi0aGxtz+vTpfPz4sUbpMjMzef36daanp7NevXr8/PPPy2N2qeTm5nLx4sXqnUltbW350UcfFRi5ysjI4N27d2lubk4AbNu2LXNzc6tNJ+/Dhw/VkxW//vprdY2uffv2/Pzzz2vUrOCqJjMzk+fOneOyZcs4cuRIenp6Ugih9vjZp08f9u3blzNmzGBQUFCZap9SfEguX76cANitW7dy9XfExcVx9OjRBMA6depw3rx5GotQamqq+pfk/v377NChQ7HLFDQlPT2d27Zt44oVK0iSvr6+fOGFF7hmzRqmpqYWmy6vNpE3NO/v709vb2/OmjWLcXFxFbKpqsjOzuaJEyc4a9YsPv/883R2dqZCoWBAQAA7duzI6dOnc9euXdXm8+gjycnJvH//Pkly/PjxbNGiBQ0NDenk5FSmfJ558YmJiaGlpSVfeeWVCneynTt3ji+//DIB0MrKiu+99x5DQ0M1bsacPn2aXl5ejIiI4JEjRzhy5EgmJCRolDY9PZ27du3im2++SRsbG7VXvdzc3HIJqkKh4NKlS/niiy/S2NiYCQkJDAoK4quvvsrly5czIiKizHnqgjzPiEuXLmW7du3UTbRmzZqRJJctW8b58+fzwIEDOhmJrClkZGSU2aneMys+UVFR/OSTT5idnc0LFy5odX1RSEgIR40aRVNTUwKgl5cX58yZw4sXL5YqRHn3V69ezQYNGjAzM5M//PADBw0axOzsbKakpKjj3Lt3j9euXWNmZibr1KlDAKxduzbHjBnDQ4cOaW34NE+8AgMD2aBBA/UUhCdPnrBfv3784osvCg3f6ytpaWk8efIkg4ODSZK9evVS92kAUHsl+OGHH7hixQoePXqUUVFR1aIfrLrxzImPQqFgYGAgbWxsaG5uzrNnz1Y4z+JITEzkihUr2KlTJ/WX28XFhcOHD+ePP/7IS5culSh6ee3nhQsXsn///oyKiqK/vz/t7OzUfTh+fn4klb/sBw4cqHSHZgqFgrdv32ZsbCxv3brFFi1aUAjB5s2bkyQ/+OADdu3ale+//z6vXbvG9PR0RkZG6sWktuJ4+PAhjx07xv/7v//jxo0bSZJNmzYtIEpvvvkmSXLixIn86KOPuGTJEu7cuVP9ay/FqeyUJD5Ceb9khBC9ASyB0p/zapILnrpvCuAXKB2GPQAwlGSk6t4MAOMA5AL4gOTBksry9fVlSEhIqTYVR1ZWFjp37oyzZ8+iS5cuWLVqFZo0aVLu/MpCXFwc9u/fj/379+PEiROIjY0FAJibm6Nly5Zo0aIFvLy80LRpUzRu3BgGBgYIDw9HaGgo+vXrh+zsbPj5+QEATE1N0aVLF4SFhaFu3bo4e/YsXnvtNdjb2+O///0vgoKCYG9vjx49eiApKQk2NjaVutdUamoqoqOj0axZM8yZMwf79u3DtWvXsH//fuTk5KBr164wMTHBwIEDsWXLFnz//fd4/Pgx6tevj1deeQWWlpZITk6Gvb09jIz0Y7s4koiOjsbNmzfxzz//oGHDhujduzd8fX1x/fp1pKenAwA++OADLFmyBPXq1YOZmRnq1q2LunXrYsiQIRg6dCjWrVsHGxsbODg4wM7ODi4uLrC1tdXxp9MPhBDnSfoWea808RFCGELpTqMngCgA5wC8TjIsX5wJAFqSfFcIMQzAIJJDhRBeADYDeB6AC4DDADxI5hZXXlnFR6FQ4MaNG9i+fTuio6Px448/4pNPPoGXlxdGjhyps83fMjIycO7cOfzxxx/4+++/cevWLTx8+BDx8fFFxm/RogX8/Pxw//59tGjRAq1atYKDgwOMjIygUCjQokULLFq0CI6Ojpg5cyaaNGmCtm3bYsuWLfDy8oK3tze2bt2K4cOHo0WLFpg5cyYWL16MRo0aYdCgQTh48CCcnZ3RqlUrREREwMbGBvb29sjJySm3GOT9gsXExGDfvn0IDw+Hs7MzPvzwQ3Tq1Al//vknAODixYu4e/cuBg4cCCEEhg4dis2bN2P8+PGIi4tDnTp1MHfuXKSnpyM4OBg2NjZo1qwZOnbsiIsXL4IkrKys0LBhQwghkJ2dDTMzs0r935JEQkIC7t27B1tbWzRu3BjTp09HdHQ0/v33X8TExGDcuHGYOnUqjI2NoVAo1GlnzpyJuXPnwsHBARYWFqhduzZsbW0xbtw4jBo1CpMmTYKFhQWsrKxgbW0Nf39/+Pn5Yffu3TA3N4eFhQXMzc3h6uoKR0dHxMXFwczMDKampjA1NYUQotI+t7apqPi0B/AlyV6q6xkAQHJ+vjgHVXFOCyGMAMRC6Uhsev64+eMVV15ZxGfatGlYvny5+heqefPmmDp1Klq3bo22bdti7dq1KOrzjRo1CpcvX8b58+ehsg8koVAo0Lp1azz//PNYvnw5cnJykJubi5ycHOTk5CA7OxsBAQE4dOgQ9uzZg/T0dDx58gSpqalITk7GyJEjMWrUKNjZ2SE5OblAmQYGBsjIyMCyZcuwbds2WFtbw9TUFAqFApmZmYiPj0dcXBwSEhKQm1usNgNQ1oosLCxgZmYGKysrPHnyBObm5nB2dsY///wDW1tbeHl54dChQ3B2dkbHjh2xfft2uLq6okePHli3bh3c3NzQp08frFy5Eu7u7ujfvz9Wr16NRo0aoU+fPvjll1/g5uaGHj16ICgoCK6urujSpQu2bdsGV1dX+Pv7Y9euXahXrx7atWuH4OBgODs7w8/PDwcPHoSTkxOee+45BAcHw83NDfXr18eOHTtgaGiIFi1aICMjAxcvXsSTJ0/w4MEDDB06FAYGBvjpp58AAH5+fvDy8sLu3bvx6NEjAMDYsWNhZGSEwMBAAMALL7wAb29vHDx4EE+ePEFOTg6GDRsGS0tLbNy4ESYmJvDz84O9vT3OnTsHIQSSkpLQt29f2NraYtu2bTAzM0ObNm1gYmKC69evw9TUFHFxcejatSvq1KmDnTt3wsLCAq1atYJCoUBERAQsLS0RHR2NDh06wM7ODjt27ICBgQHq16+PlJQUZGRkoG7dujh16hRq166trmG5urqiZ8+eWLlyJbKystSC1aFDB/Tp0wcBAQEF/s/NmzeHr68vNmzYUCC8adOm6N+/P5YvXw4AMDQ0hEKhQPv27dGtWzd88803MDQ0hJWVFUjC2dkZ3t7eCA4ORu3atWFubo64uDi4urqiW7du2Lt3L9LT02FmZoaUlBS0a9cOfn5+WLFiBYyNjeHg4ICsrCxYWVnB29sbw4YNQ5cuXUp9P4GKi88QAL1JjlddvwmgHcn388W5qooTpboOB9AOwJcAzpDcoAr/CcB+ktueKuNtKJ3Io0GDBm3v3r2r0Qd74403sHHjxkLhM2fOxLx582BsbKzeGjc/WVlZ+Pbbb/HZZ5+VOW1mZiaWLl2Kr7/+GhYWFrCwsIClpSVsbGwwcuRIjB07FtOmTYOtrS0cHBzg7OwMFxcXuLq6wtnZudRfLYVCgaSkJDx8+BAPHjxAUlISkpKSkJycjJSUFKSmpiItLQ1paWlIT09XHxkZGcjMzERmZiaysrKQnZ2NzMxMZGdnIycnB+np6WpRS09PV4tyVlaWWnglEk0YOnQogoKCNIqr9+KTn7LUfB4/foykpKRC4TY2NqhduzaKE7EGDRogOTkZjx8/VouBgYEBhBCwsrKCjY0NEhMTYWBgAENDQxgZGcHIyAjGxsY1eg/vPBEiidzc3ALXT58/feSlz3/+9N+cnBwIIdS1QCEETExMkJqaCgMDA5ibmyM5ORnGxsYwMzNDUlISjI2NYWlpiUePHsHIyAhWVlZ48OABTExMYG1tjcTERBgbG8PW1hbx8fEwNTWFra0t4uLi1OdPh5uZmRWIU6tWrQLnsbGxMDU1Re3atQudGxsbo1atWoiJiYGFhQVsbW1x7949mJqaws7ODg8fPoRCoYCNjQ1iY2NRu3ZtWFtb486dOzA1NYWDgwNSUlKQmZkJa2trxMbGwtHREdbW1ggPD4eJiQkcHR3x5MkTPHnyRB2nbt26sLKywq1bt2BqagonJyf1vvRWVlaIj49H3bp1YW1tjX/++QcmJiZwdnZGRkYGUlJS1HFcXFxgZWWFmzdvwtTUFM7OzsjMzERycrI6Tr169YqNY21tjSZNmsDBwUGj71SNbXZJJBL9piTx0eRn/ByApkKIRkIIEygdxO9+Ks5uAKNU50MAHFUNs+0GMEwIYSqEaASgKZQ+niUSyTNOqcMcJHOEEO9DuT2OIYA1JK8JIeZAOYa/G8BPANYLIW4DeAilQEEV71cAYQByAEwsaaRLIpE8O2g0z6cqkc0uiaTmUNFml0QikWgdvav5CCESAOQfprKHcgdUfUPaVTb01S5Af22rCXa5kSxyaEzvxOdphBAhxVXbdIm0q2zoq12A/tpW0+2SzS6JRKITpPhIJBKdUB3EJ1DXBhSDtKts6KtdgP7aVqPt0vs+H4lEUjOpDjUfiURSA5HiI5FIdILeio8QorcQ4qYQ4rYQYrqu7cmPECJSCHFFCBEqhNDZdGwhxBohRLzKq0BeWB0hxCEhxC3V39p6YteXQoho1TMLFUL01YFd9YUQx4QQYUKIa0KIyapwnT6zEuzS6TMTQpgJIf4WQlxS2TVbFd5ICHFW9W5uUa35LDvFuUjQ5QHlGrJwAO4ATABcAuCla7vy2RcJwF4P7OgMoA2Aq/nCvgUwXXU+HcA3emLXlwA+1vHzqgugjercGkoPnV66fmYl2KXTZwZAALBSnRsDOAvgBQC/AhimCv8RwHvlyV9faz7PA7hNMoJkFoAgAAN0bJPeQfIElAt58zMAwDrV+ToAA6vSJqBYu3QOyRiSF1TnKQCuA6gHHT+zEuzSKVSSqro0Vh0E0A1Ank+ucj8vfRWfegDu57uOgh78M/JBAL8LIc6rvDDqE04kY1TnsQCcdGnMU7wvhLisapZVeXMwP0KIhgBaQ/lrrjfP7Cm7AB0/MyGEoRAiFEA8gENQtkiSSOa5+Sz3u6mv4qPvdCTZBkAfABOFEJ11bVBRUFkv1pe5FP8F0BiAD4AYAIt0ZYgQwgrAdgAfkizgbFuXz6wIu3T+zEjmkvQB4Apli6S5tvLWV/GJBlA/37WrKkwvIBmt+hsPYAeU/xR9IU4IURcAVH+L3i6jiiEZp/oiKwCsgo6emRDCGMoXfCPJ31TBOn9mRdmlL89MZUsSgGMA2gOopfJYClTg3dRX8dHEe6JOEEJYCiGs884BvATgasmpqpT8XiVHAdilQ1vU5L3cKgZBB89MKB12/wTgOsnF+W7p9JkVZ5eun5kQwkEIUUt1bg7l9lnXoRShIapo5X9euupJ16CnvS+Uvf7hAGbp2p58drlDOfp2CcA1XdoG5Z5oMQCyoWx7jwNgB+AIgFtQ7pNWR0/sWg/gCoDLUL7sdXVgV0com1SXAYSqjr66fmYl2KXTZwagJYCLqvKvAvhcFe4OpTvk2wC2AjAtT/5yeYVEItEJ+trskkgkNRwpPhKJRCdI8ZFIJDpBio9EItEJUnwkEolOkOIjkUh0ghQfiUSiE/4fOySnCMWWou0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = [1 + i * (30 - 1) / 99 for i in range(100)]\n",
    "\n",
    "chi = pd.DataFrame({\n",
    "    'x': x,\n",
    "    'chi_1': stats.chi2.pdf(x, df=1),\n",
    "    'chi_2': stats.chi2.pdf(x, df=2),\n",
    "    'chi_5': stats.chi2.pdf(x, df=5),\n",
    "    'chi_10': stats.chi2.pdf(x, df=10),\n",
    "    'chi_20': stats.chi2.pdf(x, df=20),\n",
    "})\n",
    "fig, ax = plt.subplots(figsize=(4, 2.5))\n",
    "ax.plot(chi.x, chi.chi_1, color='black', linestyle='-', label='1')\n",
    "ax.plot(chi.x, chi.chi_2, color='black', linestyle=(0, (1, 1)), label='2')\n",
    "ax.plot(chi.x, chi.chi_5, color='black', linestyle=(0, (2, 1)), label='5')\n",
    "ax.plot(chi.x, chi.chi_10, color='black', linestyle=(0, (3, 1)), label='10')\n",
    "ax.plot(chi.x, chi.chi_20, color='black', linestyle=(0, (4, 1)), label='20')\n",
    "ax.legend(title='df')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fisher's Exact Test\n",
    "Scipy has only an implementation of Fisher's Exact test for 2x2 matrices. There is a github repository that provides a Python implementation that uses the same code as the R version. Installing this requires a Fortran compiler. \n",
    "```\n",
    "stats.fisher_exact(clicks)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stats.fisher_exact(clicks.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scientific Fraud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEYCAYAAACHjumMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAATXElEQVR4nO3de5BedX3H8feHhPsdWVIEY7AgiKNEXEALVBRBLAroKBUtZhhqOqNUvE2N1BGdqR2csSK2VpuCNF4QAUFSsQjEO7bAcpFbwCgmCoRktSAXHTDx0z/OSd2Ju8mz++z3eXaXz2tm5znnd57zO988ST77O7fnyDYRERW26HcBETFzJWAiokwCJiLKJGAiokwCJiLKzO53AZ3YfffdPW/evH6XERFjuPnmm39pe2Dj9mkRMPPmzWNoaKjfZUTEGCStGq09u0gRUSYBExFlEjARUSYBExFlEjARUSYBExFlEjARUSYBExFlEjARUSYBExFlym4VkLQ/8JURTc8BPgR8vm2fB6wETrb9cFUdEzVv0VXjev/Kc44vqiRi+iobwdi+1/Z82/OBFwO/Aa4AFgHLbO8HLGvnI2IG6tUu0tHAT22vAk4ElrTtS4CTelRDRPRYrwLmTcCX2+k5tle30w8Bc0ZbQdJCSUOShoaHh3tRY0RMsvKAkbQVcAJw6cbL3DzSYNTHGthebHvQ9uDAwB99zURETAO9GMG8GrjF9pp2fo2kPQHa17U9qCEi+qAXAXMKf9g9AlgKLGinFwBX9qCGiOiD0oCRtD1wDHD5iOZzgGMkrQBe2c5HxAxU+pWZtp8AnrFR269ozipFxAyXK3kjokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKlN4qEL2Xr/qMqSQjmIgok4CJiDIJmIgok4CJiDIJmIgok4CJiDIJmIgok4CJiDIJmIgok4CJiDIJmIgok4CJiDLVD17bRdJlku6RtFzSSyXtJulaSSva110ra4iI/qkewZwHXG37AOAgYDmwCFhmez9gWTsfETNQWcBI2hn4c+ACANtP2X4EOBFY0r5tCXBSVQ0R0V+VI5h9gGHgQkm3Sjq/fVb1HNur2/c8BMwZbWVJCyUNSRoaHh4uLDMiqlQGzGzgYOAztl8EPMFGu0O2DXi0lW0vtj1oe3BgYKCwzIioUhkw9wP3276hnb+MJnDWSNoToH1dW1hDRPRRWcDYfgj4haT926ajgbuBpcCCtm0BcGVVDRHRX9Xfyfu3wJckbQXcB5xGE2qXSDodWAWcXFxDRPRJacDYvg0YHGXR0ZXbjYipIVfyRkSZBExElEnARESZBExElEnARESZBExElEnARESZBExElEnARESZBExElEnARESZBExElEnARESZBExElEnARESZBExElKn+RruI6NK8RVeN6/0rzzm+qJLxywgmIsokYCKiTAImIsokYCKiTAImIsokYCKiTOlpakkrgceA9cA624OSdgO+AswDVgIn2364so6I6I9ejGBebnu+7Q0PYFsELLO9H7CsnY+IGagfu0gnAkva6SXASX2oISJ6oDpgDFwj6WZJC9u2ObZXt9MPAXNGW1HSQklDkoaGh4eLy4yICtW3Chxh+wFJewDXSrpn5ELbluTRVrS9GFgMMDg4OOp7ImJqKx3B2H6gfV0LXAEcCqyRtCdA+7q2soaI6J+ygJG0vaQdN0wDxwJ3AkuBBe3bFgBXVtUQEf1VuYs0B7hC0obtXGT7akk3AZdIOh1YBZxcWENE9FFZwNi+DzholPZfAUdXbTcipo5cyRsRZRIwEVEmARMRZRIwEVEmARMRZRIwEVEmARMRZRIwEVEmARMRZRIwEVEmT3aMmATjefriVHryYrWMYCKiTEcBI+kF1YVExMzT6QjmXyXdKOntknYurSgiZoyOAsb2kcBbgGcBN0u6SNIxpZVFxLTX8TEY2yuADwLvB14GfErSPZJeX1VcRExvnR6DeaGkc4HlwCuA19p+Xjt9bmF9ETGNdXqa+p+B84GzbP92Q6PtByV9sKSyiJj2Og2Y44Hf2l4PIGkLYBvbv7H9hbLqImJa6/QYzHXAtiPmt2vbIiLG1GnAbGP78Q0z7fR2NSVFxEzRacA8IengDTOSXgz8dhPvj4jo+BjMu4BLJT0ICPgT4C+rioqImaGjgLF9k6QDgP3bpntt/66TdSXNAoaAB2y/RtI+wMXAM4CbgVNtPzX+0iNiqhvPzY6HAC8EDgZOkfTWDtc7k+b6mQ0+Bpxre1/gYeD0cdQQEdNIpxfafQH4OHAETdAcAgx2sN7eNKe4z2/nRXNx3mXtW5YAJ4236IiYHjo9BjMIHGjb4+z/k8DfATu2888AHrG9rp2/H9hrtBUlLQQWAsydO3ecm42IqaDTXaQ7aQ7sdkzSa4C1tm8ed1WA7cW2B20PDgwMTKSLiOizTkcwuwN3S7oReHJDo+0TNrHO4cAJkv4C2AbYCTgP2EXS7HYUszfwwIQqj4gpr9OA+fB4O7b9AeADAJKOAt5n+y2SLgXeQHMmaQFw5Xj7jojpodPvg/kusBLYsp2+Cbhlgtt8P/AeST+hOSZzwQT7iYgprqMRjKS30Rxw3Q34U5oDs58Fju5kfdvfAb7TTt8HHDr+UiNiuun0IO87aI6pPAr//+VTe1QVFREzQ6cB8+TIq20lzQbGe8o6Ip5mOg2Y70o6C9i2/S7eS4H/rCsrImaCTgNmETAM3AH8DfANmu/njYgYU6c3O/4e+Pf2JyKiI52eRfoZoxxzsf2cSa8oImaM8dyLtME2wBtpTllHRIyp0wvtfjXi5wHbn6S5SzoiYkyd7iIdPGJ2C5oRTaejn4h4muo0JP5pxPQ6mtsGTp70aiJiRun0LNLLqwuJiJmn012k92xque1PTE45ETGTjOcs0iHA0nb+tcCNwIqKoiJiZug0YPYGDrb9GICkDwNX2f6rqsIiYvrr9FaBOcDIR4s81bZFRIyp0xHM54EbJV3Rzp9E80SAiIgxdXoW6aOS/gs4sm06zfatdWVFxEwwngevbQc8avs84P72CY0REWPq9MFrZ9N8l+4H2qYtgS9WFRURM0OnI5jXAScATwDYfpA/PEwtImJUnQbMU+1THQ0gafu6kiJipug0YC6R9G80D017G3Ad+fKpiNiMzZ5Fah9Y/xXgAJqnCuwPfMj2tcW1RcQ0t9mAsW1J37D9AqDjUJG0DfA9YOt2O5fZPrs9+3QxzUPXbgZOHfnEgoiYOTrdRbpF0iHj7PtJ4BW2DwLmA8dJegnwMeBc2/sCDwOnj7PfiJgmOg2Yw4D/kfRTSbdLukPS7ZtawY3H29kt2x8DrwAua9uX0FwVHBEz0CZ3kSTNtf1z4FUT6VzSLJrdoH2BTwM/BR6xva59y/00j6Edbd2FNI+rZe7cuRPZfET02eZGMF8DsL0K+ITtVSN/Nte57fW259PcjX0ozYHijthebHvQ9uDAwECnq0XEFLK5gNGI6Qk/osT2I8C3gZfSnOreMHLaG3hgov1GxNS2uYDxGNObJWlA0i7t9LbAMcBymqB5Q/u2BcCV4+k3IqaPzZ2mPkjSozQjmW3badp5295pE+vuCSxpj8NsAVxi++uS7gYulvQPwK3ABd39ESJiqtpkwNieNdGObd8OvGiU9vtojsdExAw3nq9riIgYlwRMRJRJwEREmQRMRJRJwEREmQRMRJRJwEREmQRMRJRJwEREmQRMRJRJwEREmQRMRJRJwEREmQRMRJRJwEREmQRMRJRJwEREmc0+2TFipHmLrur4vSvPOb6wkpgOMoKJiDIJmIgok4CJiDIJmIgok4CJiDJlASPpWZK+LeluSXdJOrNt303StZJWtK+7VtUQEf1VOYJZB7zX9oHAS4B3SDoQWAQss70fsKydj4gZqCxgbK+2fUs7/RjNc6n3Ak4ElrRvWwKcVFVDRPRXT47BSJpH8xjZG4A5tle3ix4C5oyxzkJJQ5KGhoeHe1FmREyy8oCRtAPwVeBdth8ducy2AY+2nu3FtgdtDw4MDFSXGREFSgNG0pY04fIl25e3zWsk7dku3xNYW1lDRPRP5VkkARcAy21/YsSipcCCdnoBcGVVDRHRX5U3Ox4OnArcIem2tu0s4BzgEkmnA6uAkwtriGlkPDdSQm6mnA7KAsb2DwCNsfjoqu1GxNSRK3kjokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKVH4fTLk8iD1iassIJiLKJGAiokwCJiLKJGAiokwCJiLKTOuzSBGdyhML+iMjmIgok4CJiDKVT3b8nKS1ku4c0babpGslrWhfd63afkT0X+UI5j+A4zZqWwQss70fsKydj4gZqixgbH8P+N+Nmk8ElrTTS4CTqrYfEf3X62Mwc2yvbqcfAuaM9UZJCyUNSRoaHh7uTXURMan6dpDXtgFvYvli24O2BwcGBnpYWURMll4HzBpJewK0r2t7vP2I6KFeB8xSYEE7vQC4ssfbj4geqjxN/WXgv4H9Jd0v6XTgHOAYSSuAV7bzETFDld0qYPuUMRYdXbXNiJhaciVvRJRJwEREmQRMRJRJwEREmXwfTB/ku0ni6SIjmIgok4CJiDLZRYp4GqveXc8IJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLKJGAiokwCJiLK9CVgJB0n6V5JP5G0qB81RES9ngeMpFnAp4FXAwcCp0g6sNd1RES9foxgDgV+Yvs+208BFwMn9qGOiCgm273doPQG4Djbf93OnwocZvuMjd63EFjYzu4P3NvhJnYHfjlJ5fa6/+nad3X/6bv3/Y+372fbHti4ccp+J6/txcDi8a4nacj2YEFJ5f1P176r+0/fve9/svruxy7SA8CzRszv3bZFxAzTj4C5CdhP0j6StgLeBCztQx0RUaznu0i210k6A/gmMAv4nO27JnET496tmkL9T9e+q/tP373vf1L67vlB3oh4+siVvBFRJgETEWUSMBFRZspeB9MpSQfQXAm8V9v0ALDU9vL+VdV/kg4FbPum9laM44B7bH+jYFuft/3Wye43GiPOtj5o+zpJbwb+DFgOLLb9u74WuAnT+iCvpPcDp9DcbnB/27w3zV/GxbbP6VdtnWjDcS/gBtuPj2g/zvbVXfR7Ns29XrOBa4HDgG8DxwDftP3RLvre+JICAS8HvgVg+4SJ9j3G9o6gub3kTtvXdNnXYcBy249K2hZYBBwM3A38o+1fd9H3O4ErbP+imxrH6PtLNH+X2wGPADsAlwNH0/wfXjAJ23gO8Hqaa9TWAz8GLrL9aFcd2562P+2HsOUo7VsBK4q3fVqX67+T5vaHrwErgRNHLLuly77voLkEYDvgUWCntn1b4PYu+74F+CJwFPCy9nV1O/2ySfhcbxwx/TbgNuBs4HpgUZd93wXMbqcXA58Ejmj7v7zLvn8NPAh8H3g7MDCJ/9Zub19nA2uAWe28uv37HPFv8Rrgg8APaW5G/ihN8B7VVd+T9SH04we4h+YeiI3bnw3cW7ztn3e5/h3ADu30PGAIOLOdv7XLvm8dbbqdv63LvrcA3k0zMprftt03iZ/ryNpv2vAfFdgeuKPLvpePmL5lo2Xdfi63tp/NscAFwDBwNbAA2LHLvu+k+aW5K/AYsFvbvs3IP1MX/d8xIrS2A77TTs/t9t/idD8G8y5gmaQVwIah6VxgX+CMsVbqlKTbx1oEzOmy+y3c7hbZXinpKOAySc9u++/GU5K2s/0b4MUbGiXtDPy+m45t/x44V9Kl7esaJvdY3haSdqX5zyrbw+12n5C0rsu+75R0mu0LgR9JGrQ9JOm5QLfHMdx+NtcA10jakmY39RTg48Af3Qg4DhfQ/DKdBfw9cKmk+4CX0BwemAyzaXaNtqbZBcP2z9s/x4RN62MwAJK2oNlHH3mQ9ybb6yeh7zXAq4CHN14E/ND2M7vo+1vAe2zfNqJtNvA54C22Z3XR99a2nxylfXdgT9t3TLTvUfo8Hjjc9lmT1N9KmhAU4Lbv1ZJ2AH5ge34Xfe8MnAccSXOn8ME0v5h+AbzT9o+66PtW2y8aY9mGsJ8wSc8EsP2gpF2AV9KMom/spt+27zOB04EbaD6bj9m+UNIA8FXbfz7hvqd7wFSSdAFwoe0fjLLsIttv7qLvvYF1th8aZdnhtq+faN8zkaTtgDm2fzYJfe0E7EPzW/t+22smoc/n2v5xt/30i6TnA8+jOZh+z6T1m4CJiCq50C4iyiRgIqJMAiYmhaT1km6TdJekH0l6b3sAHkmDkj7VQR8/bF/ntVerxjSXYzAxKSQ9bnuHdnoP4CLgettnT6Cvo4D32X7NpBYZPZcRTEw622tpvrD9DDWOkvR1AEkDkq5tRzrnS1rVnj5H0obbJc4BjmxHRO/uz58iJkMCJkrYvo/mwrA9Nlp0NvAt288HLqO5MHJji4Dv255v+9zaSqPSdL+SN6afI4DXAdi+WtLGFzHGDJIRTJRo785dD6ztdy3RPwmYmHTtJeafBf7Ff3wW4Xrg5PZ9x9LcwLexx4AdS4uMnkjAxGTZdsNpauA6mpv+PjLK+z4CHCvpTuCNwEM0gTLS7cD69nR3DvJOYzlNHT0laWtgvZvH17wU+Ew3NzDG1JaDvNFrc4FL2ovwnqL5UqmYoTKCiYgyOQYTEWUSMBFRJgETEWUSMBFRJgETEWX+D12bwdNtfjNrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imanishi = pd.read_csv(IMANISHI_CSV)\n",
    "imanishi.columns = [c.strip() for c in imanishi.columns]\n",
    "ax = imanishi.plot.bar(x='Digit', y=['Frequency'], legend=False,\n",
    "                      figsize=(4, 4))\n",
    "ax.set_xlabel('Digit')\n",
    "ax.set_ylabel('Frequency')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Power and Sample Size\n",
    "statsmodels has a number of methods for power calculation\n",
    "\n",
    "see e.g.: https://machinelearningmastery.com/statistical-power-and-power-analysis-in-python/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample Size: 116602.393\n"
     ]
    }
   ],
   "source": [
    "effect_size = sm.stats.proportion_effectsize(0.0121, 0.011)\n",
    "analysis = sm.stats.TTestIndPower()\n",
    "result = analysis.solve_power(effect_size=effect_size, \n",
    "                              alpha=0.05, power=0.8, alternative='larger')\n",
    "print('Sample Size: %.3f' % result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample Size: 5488.408\n"
     ]
    }
   ],
   "source": [
    "effect_size = sm.stats.proportion_effectsize(0.0165, 0.011)\n",
    "analysis = sm.stats.TTestIndPower()\n",
    "result = analysis.solve_power(effect_size=effect_size, \n",
    "                              alpha=0.05, power=0.8, alternative='larger')\n",
    "print('Sample Size: %.3f' % result)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
